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1. Introduction. The modern theories of antennas differ from the older sinusoidal 
theory in one essential point: they take into consideration the existence of the gap 
as the region from which the energy radiates, whereas the sinusoidal theory’ spreads 
the sources over the whole length of the antenna. 

The gap, in all these modern theories, is characterized by a given difference of 
potential (— V) and by this difference only. The dimensions of the gap and the elec- 
tric field in the gap do not enter the picture. If we denote by d the dimensions of the 
gap, and by £, the tangential component of the electric field along the short line repre- 
senting the length of the gap, then the physical quantity which enters the theory is: 


d/2 
-—-V= f E,dz, (1.1) 

—d/2 ; 
and the knowledge of E, and d is not assumed. Neither E, nor d appears in the thecries 
of antennas, because the gap is represented by what I shall call a 6-gap. Its model is 
obtained by the following limiting process: let us assume that the dimensions of the 
gap (d) tend to zero. Let us assume further that at the same time £, increases, but in 
such a way that the integral (1.1) remains constant. In other words, E, becomes pro- 
portional to a Dirac function: 


E, = — Vi(z), (1.2) 


where 4(z) is defined by 
6/s) = 0 for 20 (1.3) 


il 5(2)ds = 1. (1.4) 


This 6-gap model is the basis of nearly all recent investigations on the theory of 
antennas. In the final solution neither d or E, appears. They were wiped out by the 
limiting process d—0. 

Most of the investigations based on this 5-gap model employ an approximation 
procedure which is so complicated that it becomes practically impossible to carry it 
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beyond the first few steps. Even the first approximation gives a good agreement with 
experiment. But the question of convergence is hardly investigated. There is a danger 
that pushing the approximation further would spoil rather than improve such an 


agreement. 
From the theoretical point of view this is an essential point to which little atten- 


tion has been paid. 
All the existing theories of antennas can, so far as the author knows, be classified 


as follows: 

1. The sinusoidal theory which assumes a sinusoidal distribution of the current 
and ignores the gap and the boundary conditions on the surface of the antenna.! 

2. The theories based upon the model of a 6-gap, leading to an integral equation 
and its approximate? solution. In connection with these theories, that of King? may 
be mentioned, which though taking into account the boundary conditions on the 
surface, ignores the role and the structure of the gap. 

3. The theory of Stratton and Chu‘ based upon the model of a 6-gap and upon a 
rigorous solution of Maxwell’s equations. This theory is, however, restricted to an 
infinite cylinder, sphere and spheroid. The last case is the only one of any practical 
value for the construction of antennas. The advantage of the method of Stratton and 
Chu is that it avoids approximate solutions and therefore most easily allows investi- 
gation of the convergence of expressions which appear in it. 

4. The theory of Schelkunoff® which represents the antenna as consisting of two 
thin cones, and the gap as a singularity of the field at the point where the two tips 
meet. This is the only theory which takes into account the gap and does not represent 
it by a 6-model. The drawback of this theory is its great complexity and also the fact 
that it ignores the boundary conditions on the spherical ends of the cones. 

The present paper, which tries to clear up the problem of convergence and that 
of the structure of the gap, is based exclusively on Stratton and Chu’s papers, espe- 
cially those referring to the sphere and spheroid. I shall quote them as SC II and 
SC III respectively. The solutions given in SC I] and III are not based upon approxi- 
mation methods and they satisfy the boundary conditions rigorously. The restriction 
of the shape to a sphere or spheroid matters little in an investigation such as this 
which does not pretend to be of immediate practical value. It is obvious that the 
general results obtained here for the spheroid will hold also for (say) a cylindrical 
antenna, a problem much more difficult mathematically because of the boundary 


conditions at its ends. 
Any theory of an antenna leads to the calculation of 


I(0)/V = driving-point admittance, (1.5) 


where (0) is the current at the center of the gap, if the dimensions of the gap tend to 
zero. In calculating (1.5) for a sphere or spheroid, SC are led to a series 


2 E. Hallen, Nova Acta Reg. Soc. Sci. Upsaliensis 11, No. 4 (1939). 
3 L. V. King, Phil. Trans. Roy. Soc. (A) 236, 381-422 (1937). 
4 Stratton and Chu, J. of Appl. Phys. 12, 230 (1941) ; 236-240 quoted here as SC II; 241-248 quoted 
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me 
1(0)/V = Y(0) = DU Y,, (1.6) 
n=1 
by which the driving-point admittance is represented. The grave difficulty encoun- 
tered here is that the series (1.6), or rather its imaginary part, diverges. The authors 
dismiss this difficulty by a remark that this divergence is due to the 6-model of the 
gap. They write® “this difficulty is no fault of the mathematical formulation of the 
problem but results from our assumption that the voltage is applied across a segment 
of vanishing length, implying in turn an infinite field intensity and infinite current 
density at the point of application.” The authors claim further that a finite gap will 
not change the first few modes in (1.6) but will make (1.6) quickly convergent. Thus, 
according to the authors, we ought to take the first few modes for a finite gap, and 
neglect others. This sounds very vague. On the one hand if we do this, we take into 
account the fact that the gap is finite; otherwise (1.6) is meaningless. On the other 
hand, however, this finite gap will not play any role in the result. Instead of introduc- 
ing its dimensions we cross out all the modes for which is greater than a certain M. 
But how many expressions must we keep and how many must we cross out? How will 
this procedure depend on the dimensions of the gap? These are the questions which 
we shall try to clarify later. 

The same difficulty appearing in the theories of SC appears also in the theories 
listed before under 2. The difference is, so far as I know, that the difficulty which came 
out so clearly in SC, was not even discussed in the theories based on the solutions of the 
integral equations. But if the driving-point admittance (or rather its imaginary part) 
becomes infinite for a sphere and a spheroid, then it will also become infinite for a 
cylinder; this can be deduced from the integral equations of one of the theories? and 
was shown by Prof. Stevenson.’ But infinite Y(0) means zero impedance, that is zero 
resistance and zero reactance. Therefore the matching problem cannot be solved for 
a 6-gap. 

The theory gives results agreeing with experiment only because the approximate 
procedure, owing to the technical difficulties, cannot be pushed too far. It will be 
obvious from the content of this paper that any procedure which does not introduce 
the finite dimensions of the gap must lead to infinites; therefore any satisfactory 
theory must introduce these dimensions from the beginning, or, at least say in what 
stage of approximation and how they should be introduced.® 

It is the purpose of this paper to analyze the influence of a finite gap upon the 
theory of antennas. The case in which this analysis is comparatively simple is that 
of a sphere. The essential features of the difficulties just discussed come out especially 
clearly in this case. For this reason most of the paper is devoted to the spherical an- 
tenna and only at its end is the spheroidal antenna discussed and then less fully. 

This paper is in great part a result of discussions in a research group concerned 
with the problems of antennas and wave guides. In the group were Professors V. G. 
Smith, A. F. Stevenson, J. L. Synge and myself. The discussions which we have had 





6 SC III, p. 247. 

7 Privately communicated at our research group. 

® A general theory (due to J. L. Synge), which takes into account both the dimensions of the gap and 
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have been extremely helpful in clearing up many difficulties and a great part of what 
is written down here is due to my colleagues. 

To concentrate only on the logical part of the argument and not to interrupt it by 
calculations, most of them have been shifted to the Appendices. 

2. The admittance of a sphere. A metallic sphere has a radius a. A field Ey is 
applied at the surface of the sphere; we assume that EZ is a function of the colatitude 
6 and of the time ¢ through the factor e*‘. We write 


V 
Ei = — — f(0), (2.1) 
a 
where 
f {(@)dé = 1, (2.2) 
0 
and where 
V=- of Ej dé (2.3) 
0 


is the applied voltage. Our assumptions here are more general than in SC II in which 
a special f(@) was assumed. In SC II we have: 


f(0) = (= = 6), (2.4) 


the voltage being applied across a narrow strip along the equator. The admittance 
for any @ (between 0 and 7) is: 


~ = faf(0) Px (cos @) sin 0d0 
Y(6é) = ~~ = 2 sin >> Sof() (cos ) sin 


> # P, (cos 8). 2.5) 
nat ZnLf*(Px )* sin 600] ( ( 





This formula, which is basic for our further considerations is not in SC II. It isa 
generalization of SC II (19) in which f(@) =6(@) and 6=7/2 is assumed. However 
there is no difficulty in deriving (2.5) for anyone who has studied SC II. A few re- 
marks about it are added in Appendix A. 

About the notation in (2.5): P,’ are the associated spherical harmonics: 


P! = —dP,/dé. (2.6) 
The Z, are certain complex factors depending on 
ak = 2na/d (\ = wave length, a = radius) 


and the index m. We shall say more about them later, but for the moment it is sufh- 
cient to know that for very great m we have* in m.k.s. units 


Z, = — 120rin/ak. (2.7) 


* Comp. (SC II), equation (11), putting 2rwes =2m(€2/u2)“*w( eeu)? = k/60 and k =w(€2u2) 2, 
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We are interested firstly in the simple case of a gap, resembling the 6-gap, but 
finite. We assume that f(@) is a regular function, such that 


(5-9-5 +) = 


i.e. symmetric with respect to the equatorial plane. In this case only the odd n's give 
contributions to (2.5). Puttihg »=2m-+41, we can write instead of (2.5): 


° 5 £(0) sin 0PSmi1d8 
Y(0) = 2m sin >> Sof{ ane 


- : Pom41(cos 6). (2.9) 
m=0 Zama Lf* (Poms) sin 6d6 | a 





Furthermore we assume that f(@) has its maximum for @=7/2, is always positive 
and goes to zero for 8@=0 and r. For example, we could choose for f(@) the following 


function: 


2 ! 
f(0) = 2-21 @s+ 0! sin?*+! @, (2.10) 


(s!)? 
where s is an integer. Indeed for great s such a function represents a steep field, con- 
centrated symmetrically around the equatorial plane and is normalized according to 
(2.2). Though it is only a special case of f(@), it will provide us with a convenient ex- 
ample and will allow us to see what happens if s— ©, that is, if our fmite gap becomes 
a 6-gap. It is, at least, doubtful whether we are allowed to introduce immediately 
into (2.5) the 6-function and to make use of its properties. With these assumptions 
we shall calculate the input admittance for an arbitrary @ (in the interval 0, 7). This 
generalization will allow us to study the character of the divergence if s—> 0, 0-7/2. 
The calculation of /3f(0) sin 0-Psm4:d0 is done in Appendix B. The result is: 


(2s + 1)'(s + m+ 1)!(2m + 1)! 








f(0) sin 0- P2m41-d0 = 4(— 1)™(s + 1) — 2.11) 
J. ; i clea tite (2s + 2m + 3)'(s — m)!(m!)? ( 
if ms, and 
f f(0) sin 0- P2m41-d0 = 0, (2.12) 
0 
if m>s. 
Furthermore, because 
. 1 4m + 3 
| f Pom+1 sin 00 = ; ’ (2. 13) 
0 4(m + 1)(2m + 1) 


we have for the admittance (2.9) as function of 6: 
s 2s+1)!(st+m+1)!2m!(4m+3) 
Y =2n sin @ >, (—1)™(s+1) _S dK ( 
omer Zams1(2s+2m-+ 3) !(s —m) 'm!(m+1)! 


For a finite s, the sum is finite, and therefore the problem of convergence does not 
appear. Our example, as expressed by the choice of (2.10), is only one of the possible 
examples. We could have taken a properly normalized step function instead, as we 





Pomsi(cos 0) (2.14) 
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shall do in Appendix E. Our choice (2.10), however, has the advantage that the series 
by which Y is represented is finite. 
Let us write 


Y = Y(6,s), (2.15) 
indicating the explicit dependence of Y on @ and s. First we wish to calculate: 
Y(6) = lim Y(@, s). (2.16) 


To do this, we split Y(@, s) into two parts: 
Y(6, s) = Vo,m(0, s) + Yar.(9, 5), (2.17) 


where: 


M 2s+1)! 1) !2m!(4m+3) Pom. 6 
tiie tai) eons) 
Zom+1(25-+2m-+ 3) !(s—m) !m!(m+-1)! 





m=0 


Yu,s(0,s)=2rsin@ >> as above. (2.19) 


m=M+1 


Here M is assumed to be an arbitrary but great number, later kept constant, 
while so The Yar,.(9, s) is calculated approximately by the use of Stirling’s formula 
for m!, and by replacing summation by integration. The calculations are performed 
in Appendix C and the result is 


iak\/sind ~ 7 3 
Yu,.(9, 5) = —————— > (-—1)"' m™ cos ES + (2m 2 5) 


607 m=M+1 
m/s 1+ c 
“exp | - sf log —— ar|. (2.20) 
0 i-f 


Now if s—o, we have (App. C), 


“ 3 
> (—1)"# mr cos E + (2m + ~)o]. (2.21) 


T m=M+1 ~ 


iaky/sin 0 
Yu.) = — nes 





Let us now investigate this expression near the equatorial plane by putting 


rT d 7 d 
6=— -- or @=—+—,» (2.22) 
2 2a 2 2a 


where d is small compared to a. Then we have: 


d iak ¢™ 3\d 
V0 (-) ~ — f m~* cos | (2m + =) <| dm, (2.23) 
2a 607 J wv 2/2a 


and, for small d, we finally obtain the leading term: 


d tak =—6f( Md 
rai) =- tots 
2a 607 a 
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We see the character of the divergence of our series if s—«, d—-0: Yy,..(—0) in- 
creases to infinity like log d/a: 
| Veow| ~ | log (d/a) | . (2.25) 


Therefore we cannot define the driving-point admittance by a 6-gap with @=7/2. 
Such a driving-point admittance does not exist. There are three possible ways of sav- 
ing this situation, which we can represent schematically in the following way: 

1. f(0) +5; d¥0. 

2. f(0) +5; d=0. 

3. f(0) =6; d¥0. 

In the first case we take a continuous field distribution and Y at the distance d/2 
from the equatorial plane. The difficulty in accepting this view is that we must know 
f(@) and connect d with f(@) in some way. The most obvious way would be to assume 
that f(@) is a step function, and to identify d with the width of this step function on 
the one hand, and with the dimensions of the gap on the other. This special choice of 
f(@) is discussed in Appendix E. Generally, however, we must introduce three arbi- 
trary factors: (a) the function f(@), (b) the dimension of the gap, (c) some connection 
between d and f(@). As long as there is no experimental evidence to guide us in deter- 
mining all these factors, this view seems to be too arbitrary for a satisfactory theory. 

In the second case the driving-point admittance is given by (2.9) for @=2/2. Here 
only the function f(@) is arbitrary and on its choice Y will depend. 

In the third case we go back to the 6-function but calculate Y, not in the equatorial 
plane but at the distance d from it. In this case the definition of the driving-point 
admittance will be simplest if we can identify d with the physical dimensions of the 
gap. We shall accept here this point of view, assuming that d—the only new parameter 
entering the theory—represents the known dimension of the gap. 

It is possible that one of the other views will turn out to be more convenient and 
the present discussion tends to show how the frame of the theory of the antenna 
must be broadened to include the size of the gap. 

At present however, we shall accept the simplest view: 6 electric field and finite gap, 
corresponding to our third case. 

For this case Yy,.(d/2a) is expressed in (2.24) and the remaining task is to cal- 


d d 
Vo.m (s. -) = vou(=). (2.26) 
ee 2a 2a 


This is done in Appendix C. According to (C.2), (C.3) and (C.9), we have: 


1 1)!(2 1)! 1 
(Gt I)(s + m + 1)!(2s + . (2.27) 


culate 


s7@ (s — m) !(2s + 2m + 3)! hel Q2m+2 


Therefore 


d 2m!(4m + 3 > d 
Y ( ) -- or cos (< )> > (-— 1) net — ) — Poms 1 (sin <) 
2a 2a ” Deme2 2(m a. 2a 


m=0 


iak Md 
aio c (2.28) 
607 a 
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which we interpret as the final result representing the input admittance at the end 
of a small but finite gap in the center of a conducting sphere, when the field in the 
gap is represented by a 6-function. 
For d—0 we have, because 
(2m + 1)! 


Peni1(0) = (— 1)*»—_——__ 2.29 
om+1(0) = (— 1) (m*2°™ ( ) 


the following formula for Yo,u(d->0) 


M — [P2m41(0) ]?(4m + 3) 


T 
Yo.w(0) = — } 


(2.30) 
2 m=0 Eins i(m aa 1)(2m “ 3) 





This is exactly the formula for the driving-point admittance obtained by SC II, (19), 
if we there replace nm by 2m+1. Though our discussion shows that Yy,..(0) and there- 
fore Y(0) becomes infinite at the driving-point, one could argue, as SC does: 

In Yo,u(d/2a) the gap d will play a small role and therefore we may assume d=0. 
On the other hand, if d is finite the integral representing Yy,,..(d/2a) can be neglected 
and thus we arrive at a formula (2.30) which is convergent and does not contain the 
dimensions of the gap. But the argument could just as well be reversed. We could say: 
For small d the end of the series represented by the Ci function becomes essential, 
whereas a large d will modify even the first few expressions of our series. In a satis- 
factory theory we must be able to answer the following question: how can we choose 
M for a given ak, d, €>0 so that 


| Yar. | <e. (2.31) 


Our discussions provide an answer to this question which was not supplied by SC. 
Let us take as an example: 
ka=v/2; d/a=0.1; M = 15. 


Then the absolute value of the first term in the sum in (2.28) is approximately 
3n ait 
—|Z; |. 
2 

But, as will be seen from equations (3.5) and (3.16), 


—1 | /2 


2x X 60 | 
therefore the first expression is 34/2/4 X60. On the other hand the absolute value of 
the integral expression is (1/2/60) Ci(1.5)~(4/2/2m X60), which is of the same order 
as the first expression in the series (2.28). 

The result as expressed by (2.28) could have been found much more quickly by 
the use of Dirac’s 6-function. Indeed, putting f(@) = 5[(1/2)—6@] in (2.9), we find 





d 1\ 2. Pomy1(0)P2m41(sin d/2a) 
r(=) = 29 cos (~) Fe sor ROD Fomish (2.32) 
2a 2a] mut Zom+t [/*(Pimss)” sin 6d6 | 


By splitting this into two parts and using the asymptotic formula for P2m41, we 
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again obtain (2.28). But in obtaining (2.28) by the transition s—«, we avoided the 
use of Dirac’s 6-function and at the same time gave an example of the input admit- 
tance with a given field distribution in the gap. Though the choice of M is arbitrary, 
as long as we can apply Stirling’s formula, our equation (2.28) will be the more exact 
the greater is the chosen M. 

3. The coefficients Z:2,,,;. The coefficients Z, are defined in SC II, (8), (6). Assum- 
ing ¢ (conductivity) —«, we have, in m.k.s. units, 


(2) 
H,~1/2( ka n 
Z.= i 20-60) - val  — 


i te (3.1 
He 2(ka) = ka , 


We shall now show how these expressions can be calculated without the knowledge 
of Hankel’s functions. The reason for doing it is that we intend to apply a similar 
method later when dealing with the problem of a spheroidal antenna. 

Because of the recurrence formula: 





9? 
(2) 2m+1 _ ¢@) (2), 1 (2) 
Homi1/2 — — Hom+sj2 = Hom+sj2 + —— Hom+s/2 (3.2) 
ka 2ak 
we can write (3.1) 
4 (x !2Homasy2)’ 
Lom+1 _ t- 27-60 Es 7 ae (3.3) 
n/2 m+3/2 
where we put 
x = ka. (3.4) 


But the x'/2HY.,. are “densities” corresponding to Hankel’s functions of the second 
kind. We write 
(2) 
x? om 32 = hom+3/2 => h, (3.5) 
where hom 3/2 satisfies the equation: 
4m? + 6m + 2 
h” + : h+h=0. (3.6) 


x 





It is more convenient to calculate Zen4: directly, and not by means of the h’s. Indeed 
we can easily find the differential equation which Z must satisfy. We write: 





hom 3/2 
Somp1 = — - (3.7) 
hom+s/2 
2x X 60 
Zot = §-—— (3-8) 
Som+1 
Then (omitting the indices) 
bed (=) h”’ (3.9) 
.s wv) h- 


which, because of (3.6) and (3.7), becomes 
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4 (= - 1) oe (3. 10) 


x 
where 
vy? = 4m? + 6m + 1. (3.11) 
Riccati’s equation (3.10) can be solved by a recurrence formula. Its derivation is 
given in Appendix D. This.formula is: 
Com—1 X [x?/(4m + 1) — 2m] + x? 


ee Sette elt, SS on 1 -» (3.12) 
Som—1[2m(2m + 1) — x?] + x} [x?/(4m + 1)] — (2m + 1)} 


and it is a purely algebraical formula which does contain derivatives. We start by 


finding ¢_,. For m= —1, equation (3.10) takes the form: 
’ 2 > 
- wie rs = 1, (343) 


As a solution of (3.13) we choose 
fi =i. (3.14) 


In justifying (3.14) we must go back to (3.5) and (3.7). If 
h’ 
—~-—-i for x7 © 
h 

then 


e2A?) wh~e-* for x7 @, 


We can now calculate the first ¢ in our series corresponding to m=0. We have, from 
(3.12) and (3.14), 


——; (3.15) 


or, because of (3.8): 


1 1 x ix ; 
i a a. | a 
Zi 2r X 60 L1 — x? + x4 1 — x? + x4 


We see: for x—>0, the real part of {; goes to zero as x and the imaginary part as x‘. 
For x, the real part of ¢ goes to zero as 1/x* and the imaginary part goes to 1. 
From the recurrence formula it follows that for all m and small x(x<2m) the real part 
of ¢ goes to zero as x, and is negative, whereas the imaginary part goes to zero at 
least as x4. 

For very great m and very small x, (x2m) we find the leading term by writing 
(3.10): 


. 4m - 
i+ °—= y (3.17) 


x? 


The leading term is: 
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Sompr = — 2/2m> (3.18 


in agreement with (2.7), if we take into account (3.8). 

4. The spheroidal antenna. We shall now generalize our arguments for the case 
of a spheroid. 

We characterize the points in space by the spheroidal coordinates £, 7, ® with: 


&2 1; —-13839751; OS @8 2, (4.1) 


where &, , ® are defined in terms of cylindrical coordinates by 
P= P(e —1)\(1—7); 2 = fen; $= ¢, (4.2) 
and in terms of polar coordinates by 


én 
(B+ 9? = 1)" 


Here 2f is the distance between the foci.* A spheroid is characterized by a constant 
coordinate §=£. Knowing £ and f we can find the half-axes a, b: 


a= fo; b= fVP—T. (4.4) 


= 9, (4.3) 





R= f(%+7—1); cosd= 


In CS II], the input admittance of such a spheroidal antenna is calculated. Again, 
as in the case of a spheroidal antenna, the 6-model of the gap is assumed and the 
imaginary part of the input admittance is infinite. 

We start by generalizing CS III in two respects: first, we do not take the 6-model 
of a gap but an arbitrary field; secondly we write down the admittance not for 7»=0, 
but for any 7 (|| $1). The argument is very similar to that in the case of a sphere 
and hardly worth repeating. The result of these generalizations leads us to the follow- 
ing formula corresponding to (2.5) in the case of a sphere: 


© Til — 2) fln)Ser 
y = 21-2) 5 S-a(1 = 0°) fl) Ser(n)dn 


Sei(n). 4.5 
5 ZifeCt = wSe(a) Pan (4.5) 





It is written in such a way as to expose an analogy with (2.5) and, at the same time, 
to use the notation in SC III. To establish a connection with SC III one must sub- 
stitute in (4.5) 

; 2 4 
Qn ” tkf (to sents 1) Re;(o) (4 6) 
Zi 60 d/dto[(e — 1) Re*(é)] . 





Equations (4.5) and (4.6) are, for 7 =0 and f(n) = 6(n), identical with SC III (2.8). 
For the sake of completeness we shall explain the notation in (4.5) and (4.6): 
The function f(n) satisfies the condition 


f(n)dn = 1, (4.7) 


* This f will hardly be confused with f(n) characterizing the electric field and always written with its 


argument 7. 
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and in the case of a symmetrical field: 


f(n) = f(— »). (4.8) 
The functions Se} and Ref are solutions of the equation: 
(1 — 22)W” — 42W’ + (A — f7h*2?)W = 0 (4.9) 


for =n and z= respectively. The eigenvalues of \; are determined by the condition 
that Se remains finite at the poles n= +1. The function Re for given eigenvalues 
of \; is so chosen that the field at large distances from the center reduces to a wave 
travelling radially outward. The functions Se} and Re are related to the spherical 
harmonics and Hankel’s functions through coefficients which must be determined 
numerically. However, for our argument the knowledge of these relations is not es- 
sential. 

The prolate spheroid has a gap across the equator 7=0. By a method similar to 
that of Section 2, we shall now calculate the input admittance at the points 


s= + 4/2, (4.10) 
which, because of the second equation in (4.2) means: 
n = + d/2fto, (4.11) 
or, finally, because of (4.4): 
7 = + d/2a, (4.12) 


where a is here half the long axis. For f(y) we shall take the 6-function: 
f(n) = 5(n), (4.13) 


and our problem will be to find what happens to (4.5) if d—0. It would have been 
more rigorous to take, say: 

f(a) ~ (1 — 9°)° (4.14) 
and then go to the limit s— «, as we did in Section 2. However it will be much quicker 
to use the 6-functions instead and we have seen in Section 2 that the result is the same, 
which is also true in the case of a spheroid. 

We now rewrite (4.5), introducing f(y) = 6 and 7 =d/2a. Because Se; =0 if 1 is odd, 
we have: 
1 1 
S€2m(0)Se2m(d/2a) 


v(=) = 2x[1 — (d/2a)?] >> —— ; (4.15) 
2a 


m=0 ZomJt\(1 - n*)(Set_)*dn 





Now again we write (4.15) in the form 
VY = You + Yun; (4.16) 
and shall proceed to calculate Yy,,. approximately for sufficiently great M. For the 
moment we shall ignore Z2,, and shall calculate all other expressions appearing in 


(4.15). 
The functions Se,, must be finite and satisfy (4.9) for z=7. But if Xen is very great, 
which we assume, f?k?y can be neglected in (4.9), because | n| <1. Therefore, for great 


om We can write (4.9): 
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a= n’)Seom — AnSeom + NomS€2m = 0. (4.17) 

Therefore, Se}, will not depend on f for large Ae. Substituting in (4.17): 
Seim = (1 — 0°)-¥?P, (4.18) 


we can write instead of (4.17): 


9 





P + AemP = 0. (4.19) 


2 


n 
P” + 2nP’ — 
.=¢ 


But this is exactly the equation for Ps'm4: for which an everywhere finite solution 
exists if 


hom = 2m(Qm + 3) ~ 4m?, (4.20) 


and m is very large. The approximate expression for Pm41 with large m is: 


gt BON 8 = I(2 4 ~) =| (4.21) 
tnaa(5) = 2 "WY cemltlia WN” BABS 


therefore, because of (4.18), 


1 1 d d\? 
S€2m(0)Seom (—)[ (-) | 
2a 2a 
wood FE 1 (2 me ~) =| (4.22) 
*f-QTeale+ Dg « 


The next step is to calculate the integral in the denominator of (4.15). Again be- 
cause of (4.18), we have: ; 


on ae 4(m + 1)(2m + 1) 
f (1 ae n°) (Seom)2dn = f (Pom+i)*dn = ” M ~~ 2m. (4. 23) 
1 a | 4m + 3 





We can, therefore, write for small d and sufficiently great M: 


— [2 = cos [(2m + Be reN 





(4.24) 


p> 


T m=—M+1 Zom 
The last step will be to find Zom. 

5. The coefficients Z2,, for a spheroid. It will be a little more convenient to intro- 
duce, instead of Re, the corresponding density-functions, as we did in Section 3. 
The equation which Ré&,, satisfies is: 

” , . 4 
(¢? — 1)Ress’ + 4¢Resn, — (Nom — {2k Reom = 0. (5.1) 
We now introduce 
(2 — 1)Ream = P, (5.2) 
and (5.1) becomes, because of (4.20), 


(—? — 1)P” — [(4m? + 6m + 2) — frke?]P = 0. (5.3) 
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Because of (4.6) we have: 





= =- “ P/P’ = — — Em (5.4) 
where 
P/P’ = fom = £. 
We have: 
ti = 1 — g2P”/P. (5.5) 
Putting 
4m? + 6m +2 = <> pk? = c’, (5.6) 
we get from (5.5) because of (5.3) 
Cf ee en ne 
e#-1 @-1 
which is Riccati’s equation for ¢. Introducing now: 
p= p® — ¢?, (5.8) 


we can write (5.7): 


Pp" 
y+ “(5 a “*) = 1. (5.9) 
ie 


The leading term of the solution for every great 9 and small & (so that p?/&? —1>>c?) 


P 1 as 1 os 
— = fom ~ + —JVBP-—-1~\~+—V#—- 1. (5.10) 
p’ D 2m 


fre} 


The proper sign will be selected by the condition that for £0 and f—0, our solu- 
tion goes over into that of a spherical antenna. 
Introducing (5.10) into (5.4) and into (4.24) we have: 


d 2ikf 2. VB-1 [ ( 3\d 
Trait) = + — — cos | { 2m + =) | 5.11) 
r (<.) 607 m 2 2m 2 2a ( 


Changing the summation into an integration we have, because of (4.4) 


d kb ¢* cos |(2 3) (d/2 
ua(a)- taf SE (5.12) 





607 M ™m 


But for b=a this is exactly (2.23) if we take the lower sign in (5.10). Therefore finally: 


d ibk (Md 
tela Sl 19) 
2a 607 a 
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and because of (4.15): 


1 1 
d d\?] & S€2m(0)Se2m(d/2a ibk Md 
r( ) = an 1 -- (<) | z onl (4/24) — ci(=*). (5.14) 
2a Saf Jano Zomf (1 —n°)(Se! )%dn 600 a 
As we see, the only difference between (5.13) and (2.24) is that in the case of the 
spheroid we have as a factor the small axis. Thus Yuy,,..(d-—>0) tends to infinity like 


b log (d/a) (5.15) 


and can be made to be always finite if b log d/a is finite. 

For thin antennas and finite d, Yy,,. is much smaller than in the case of a sphere 
with radius a; therefore the neglect of Yy,,.. is much more justified. This is however 
not true for thick antennas. But the essential result is the same as in the case of a 
spherical antenna: the dimensions of the gap must be involved in the theory of the 
antenna, otherwise the input admittance becomes infinite. But infinite driving-point 
admittance means zero driving-point impedance. In view of this result it is difficult to 
see what the curves 3 and 4 in SC III mean. They seem to represent the resistance 
and reactance of a spheroidal antenna and a current distribution which does not be- 
come infinite for 7 =0. They may refer to a finite gap and to an electric field in the 
gap, but in that case we would like to know to what field and to what gap they refer. 
[t is plausible that neither the field nor the size of the gap matters much if the antenna 
is thin and if the size of the gap is greater than the thickness of the antenna and much 
smaller than its length. But such a statement ought to be deduced from the theory or 
at least formulated explicitly as an assumption. 





APPENDIX A. ON THE DERIVATION OF (2.5) 


Equations (6) and (8) in SC IT give 


co be 
Hejraa = — >, — e**'P,! (cos 8), (A.1) 


n=1 n 
where A, is here used instead of (wu:/io;)'/2C, in SC’s notation. On the other hand 
SO Pa CF 


) is: 
V 7. 
— — f(0) = Ef = et)” A, P,! (cos 8). (A. 2) 
a n=1 


From (A.2) it follows because of the orthogonal properties of P,’ , that: 


ek De lf (0)(P,!) si oi ownd (A.3) 
—— 2n(n + 1) 0 J ecmet a ; 


Putting (A.3) into (A.1) we have: 


ee 2n + 1 r P! 
Ons oes ee f 0)(Px) sin oio | 
es a 2n(n+ =| i f(0)(Pr) ( a 


V = Sof(0)(P.t) sin odo 


2 


a wat Znf*(Px)? sin 640 





P,! (cos 6). (A.4) 
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Finally from SC II (13) we see: 


= 27a sin 0H e\R=a} 


therefore 
tf *  (5(0)(P,!) sin 6d@ 
Y = — = 2rsiné>, Jf O(Pr) : 
V ZnJ*(Px )? sin 6d0 


n=1 





P, (cos 0) (A.5) 


which is exactly (2.5). 


APPENDIX B. THE CALCULATION OF 











Qs+ Di f°. ’ 
B = 2-2:-1 — f sin?*t? §Pom41d0. (B.1) 
(s!)? 0 
Because of 
| ae = dP om+1/d0 (B.2) 
and the recursion formula 

2(m + 1) 2m + 1 
COS OP amp. = ———— Pam) + Pom; (B.3) 

4m + 3 4m + 3 


we have, introducing uy =cos 6: 


td +1 
f Phmsi sin**+? 9d0 = 2(s + »f Pawis(1 — 2*)*udp 
0 1 


2s +1) +1 
= ‘mes 420 + vf (1 abi B?)* Po ¢m+1du 
+1 
+ (2m + 1) (i - 1*)"Poad\ (B.4) 
anil 


Therefore, the problem is reduced to finding an integral of the type 
; +1 
T= f (1-19) "Pad (B.5) 
—1 


To find T, we must combine the results of Examples 1 and 6 on pp. 310, 311, in Whit- 
taker and Watson’s Modern Analysis, Fourth Edition, 1940. The result is: 





‘ !\? 2n)! ! 
T, =(- raze (=) - (2n) = n) Ss te, 
n!J} (s — n)!(2s + 2n + 1)! (B.6) 
T.=0 if s<n. 
The combination of (B.1) and (B.4), (B.5), (B.6) gives: 
s+ 1)(2s + 1)! . . 
B= 2-* ( \ [2(m + 1) Timor + (2m + 1)T mI, (B.7) 





(4m + 3)(s!)? 
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which, because of (B.6) and some simple algebra, gives the results as expressed in 
(2.11) and (2.12). 


APPENDIX C. CALCULATION OF Yy,,,(9, s) 


Because of (2.7) we can write (2.19) 





ee 
peter Ca 
2» 1 iak sin 0 
i (25-+-1) (s-+m-+1) !2m!(4m+3) 
aa et ae P2m @). (C.1 
' 2 | S41) Gs am43)s—m) mimi) m4) ae. &H 


First, we calculate the expressions in which s appears. These are: 
(s+1)[1-2-3--- (s+m+1)][1-2--- (2s+1)] 
[1-2-- + (s—m)][1-2--- (2s+1)(2s+2) --- (2s+2m+3)] 
[(st+1—m) --- (s+1+m—m) ][(s+1)(s+2) - -- (s+m+1)] 
22m+2(s-4-1)(s+1+4+4)(s+2)(s +244) - ~~ (s+-m+1)(s+m+144) 
[1—a/(s+1)] 


1 m 
phinrerone MN ghecceteet eetecs 4 (C.2) 
geme? 25 [1+(at+4)/(s+1) 











Assuming that s is great, taking log of the above expression and changing the 
summation into integration, we have: 


1 m 1—a/(s+1 1 = 1 
II We A 1 - exp | - sf log ; +é 
0 


22 ot [1 4 (a + 1) /(s + 1)] om ‘Q2m+2 i 








ie], (C.3) 


where 
x=m/s- (C.4) 


If for the expressions depending on m, we use Stirling’s formula, we obtain: 





r 1 * 1+é . 
C,, = ——= (— 1)" exp |- sf log dé | Poms1(cos 8). (C.5) 
2mv/ rm 0 i-—§€ 


Now for Pe'm41 with great m we have approximately: 


, m T 3 
Pomyi = — 24/ cos (= + (2m oe -)*), (C.6) 
3 sin 6 4 2 
therefore 
° (— 1)™*! f° i+é | [- ( ~)e| 
C. = —— —s§ lo dé |}cos| —+ | 2m + — 6], (C.7) 
mrv/ sin 6 ai | 0 . 1-é€& : 4 2 


and because of (C.1) we have (2.20): 
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takv/ sin 6 2 1 
Yu,s =——_ © (- 1)" — 
607 m=M+1 





=f, pale +5) 0 

exp| — s ) cos | — m+ — ‘ ; 

. 0 . i-¢ 4 2 

If s—> the expression in the exponent is: 
m/s 1 oe é 1 ma 1 -f- 

lim | - sf log as] = — lim |—f log at = 0. 

se 0 1 = E a0 Qa 0 


hel 
Therefore (2.21) follows from (C.8) and (C.9). 








(C.9) 


APPENDIX D. THE RECURSION FORMULA (3.12) 


We can deduce the recursion formula (3.12) for any cylindrical function density z, 
satisfying the equation: 




















as 
sf? — — 2p +2, =0 (D. 1) 
The known recursion formula for 2, is: 
* eae 3 ’ 
: Zp-1 — Zp-1 = Zp. (D.2) 
x 
On the other hand we have, because of (D.2), 
Lie a ae (p — 4)a-%y1 + (p — 4a tp — 25-1 (D.3) 
Zp (p — 3) Sp — 25-1 
This gives, because of (D.1): 
Sp — (p — ¥en a + (p — $a sn  — {[(p — 1)? — d)a? — 1} 2,4 :4) 
Zp (p — $)x'*sy1 — 55-1 ; 
Putting 
tp /Zp = Yn (D.5) 
we have 
1+ (p — Baye a — (p — DPE 
Yp = —— (D.6) 
(p — $)a7? — yo 
Now, by employing this formula twice, we have: 


7 — 


Yp—2 x(x?/(2p — 2) — (p — 4)) + (p — 3/2)(p — 4) — x? 
Yp—-2x* + x(x*/(2p — 2) — (p — 3/2)) 
Let us now introduce 





p = 2m + 3/2, (D.7) 
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as dedicated by the comparison of 27 /z, with (3.4). Thus we have: 


x (x?/(4m + 1) “— (2m + 1)) yom—1/2 + 2m(2m + 1) — x? 
V2m+3/2 = — . (D.8) 
x Vom—1y2 + 2%(x?/(4m + 1) — 2m) 





Finally, introducing: 


1 
’ (D.9) 





Comp = 
V2m+3/2 


we have our recursion formula (3.12). 


APPENDIX E. THE DRIVING-POINT ADMITTANCE IF f(@) IS A STEP FUNCTION 


We shall calculate the driving-point admittance for a spherical antenna, assuming 
for f(@) a step function defined by: 


(0) : f . + f 262 vd . 
es. Ss ae ee ae 
(E.1) 
ish: thie ib ap ee aay aa hai cae 
= or —+— an ———-: 
f 2 2a 2 2a 
We assume that d is small and we shall calculate the admittance for 
d T d , 
cos #6 = —» thatis #=—--—- (E.2) 
2a 2 2a 


We introduce (E.1) into (2.9). We have, because of (2.6) and because d is small: 


, d\ a 
f S(O) Poms sin 6d6 = 2Pom+t (~)< : (E.3) 
0 2a/ d 


Therefore (2.9) can be written for our choice of f(@): 








(2 4a. Pomy1(d/2a)P2m41(d/2c) (E.4) 
2a) dd. Amo Lf ( Poms)? Sin 040 |Zemss 
We are interested in calculating Yy,.(d/2a) defined by: 
d\ 4 © Pamsi(d/2a) Pom41(d/2a 
Pun(—) cick 2 sal! ®) 2 +a(d/ ) ’ (E.5) 
2a d m=M+1 Zom+1|fo(P2m+1)* sin 6d6 | 


For great m we have: 


d Ea 3\d 
Panss( =) = (— n4/— sin [ (2m + =] 
, d m 3\d 
Plans(33) = (- 1n24/~ co [(2" * at 
| f " (Phnes)? sin oi eae 

0 2m 


(E.6) 
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1 tak 
Zemsr 4 X 6Om 

















Therefore 
; d iak = sin [(2m + 4)d/a] 
Vit, (~) OD pcemeres a - : (E.7) 
2a 60 ae w M+1 2m*d/a 
Changing summation into integration we have: 
d iak = sin x 
Yua(—-) = —— f dx, (E.8) 
2a TX 60 2M (d/a) x? P 
and the leading part of (E.8) is for sufficiently small d: 
; d iak _{2Md 
Vt. (=) =— - Ci : (E.9) 
2a 3 X 60 a 


It is interesting to compare this result with (2.24). The only difference is in the 
argument of the Cz function, which matters very little, if the argument is small. 
One would be tempted to calculate the driving-point admittance by taking 0 =7/2 
and not, as we did, @=2/2—d/2a. This means that in (E.4) and (E.5) we would have 
to substitute zero instead of d/2a as an argument in P?n41. Then for Yy,..(0) we easily 
find: 
iak = sin [(2m + 3/2)d/2a] 


Yu..(0) =- 
ml) x X 60 p> m*d/a 








(E. 10) 
and changing summation into integration we have 


tak Md 
Yy..(0) = — ci(=*), 
r X 60 a 





that is the same result as in (2.24). . 
The argument will be very similar in the case of a spheroid. Thus we can collect 


all the results concerning Yy,,.. in the following two tables: 












































Sphere Spheroid 
— = ae —_—- a ae ee a a 
| | | | 
= 00 | 0 00 
Dirac’s 2 ee Dirac’s apes 
function | : || function d 
| 27 2a iak (Md || | +36 | ibk . (Md 
A -o(t) ———# | a(t 
hl | 607 a | 0 60x a 
Step a 2 | Step 
function x d iak _. (2Md\ || function d ibk _. (2Md 
—+— —— (—) ] +— _— ci ( ) 
2 2a 607 a i j 2a 607 a 
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MULTIPLE REFLECTIONS BY PLANE MIRRORS* 


BY 


L. B. TUCKERMAN 
National Bureau of Standards 


The article by Synge’ called to my mind a more general treatment of reflections 
by plane mirrors which I saw when, in 1898, as a sophomore in college, I was studying 
quaternions. I cannot cite the reference because I took no particular note of it at the 
time. I never had occasion to use the theory until 1921 when I needed it to determine 
the errors in my optical lever system,? which consists of an autocollimated three- 
mirror combination. 

My attention has since been called to a book by Silberstein* which treats the same 
problem with a modified Gibbs notation. Silberstein’s Eq. (25) on p. 22 is identical 
in form with Synge’s Eq. (4.1) but his later treatment is different and he does not 
note, a Synge does, that the formulae represent rotations. 

The quaternion treatment of the problem is so simple that I had no trouble in 
reconstructing it, at least in its essential features. Comparing it with Synge’s and 
Silberstein’s treatments strengthens my feeling that Heaviside, Féppl, Gibbs and 
their followers did the world a disservice in rejecting Hamilton’s simple, concise nota- 
tion and algorithms for the cumbersome complexity of the vector notations currently 
taught. It therefore occurred to me that the readers of the Quarterly of Applied 
Mathematics might be interested in seeing an outline of the quaternion treatment. 

When the first draft of this paper was submitted for publication, the reviewer 
called attention to a recent article by Coxeter.‘ He gives a very general quaternion 
treatment of reflection and rotation in three and four dimensions but does not dis- 
cuss special cases in detail. 

In his introduction Coxeter states: “Apparently none of these men thought of con- 
sidering first the simpler operation of reflection and deducing a rotation as the product 
of two reflections.” 

The article which I read, presented a rotation as the simpler operation and consid- 
ered a reflection as a rotatory inversion, compounded of a rotation and a central 
inversion. This is the order used in the following discussion. 

The underlying basic theorem of reflection at plane mirrors is not confined to rays 
of light incident upon the mirror but applies to the whole image space. It may be 
stated as follows: The image space formed by a plane mirror can be derived from the 
object space by rotating the object space 180° about any normal to the plane, followed by 
a central inversion about the foot of the normal. 

In the quaternion notation the square of a vector is negative instead of positive. 
This, as Hamilton showed, makes it possible to treat a “quaternion,” the sum of a 
scalar and a vector, as a single quantity in addition, multiplication, division and 
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differentiation instead of introducing the multiplicity of scalar and vector operations 
of the notations of Heaviside, Féppl and Gibbs. 

Computations with quaternions are thus carried out with the simplicity of ordi- 
nary algebraic computations. The only complexity, if it may be called a complexity, 
is the necessity, common to all vector notations, of remembering that in general 
abba. The separation of the quaternion results into their scalar and vector parts 
can be made at any desired stages in the process by the operators, S=scalar part and 
V =vector part. 

The quaternion formula for rotation of a vector i through an angle # about a unit 
vector n as axis is: 


i’ = gig”, (1) 
where the quaternion 
q = cos 40 + n sin 30 (2) 


is the sum of Sg=cos $0 and Vg=nsin 4. 
If # =180°, this reduces to: 


i’ = nin. (3) 


The central inversion merely requires changing the sign of all vectors passing 
through the foot of the normal, giving for the effect of a single plane mirror 


i’ = — nin“ = nin. (4) 


Equations (1) and (2) are identical in form with Coxeter’s (p. 139) Theorem 3.2 
and Eq. (4) with his Theorem 3.1. 

The equivalence of Eq. (4) with Synge’s Eq. (4.1) may readily be verified by ex- 
pressing it in quaternion notation. Writing (n-i) = Sni and n?= —1, gives: 


ii =i — 2n(m-i) = i+ 2nSni =i+n(ni+in) = (1+ n%i+nin=nin. (5) 
The effect of two mirrors is obtained by simple iteration: 
i” = nei’n, = neniinym, = gig’. (6) 


The second central inversion cancels the first, so the result is a pure rotation. The 
angle # and the direction of the axis a are given by 


cos 308 + asin 30 = g = nom; = Sno, + Van. (7) 


Here Snon; = —cos (me, ni), where (me, ni) represents the angle between the normals 
to the two mirrors, and Vnen is a vector with a length equal to sin (mz, mi) per- 
pendicular to nz and n;, and therefore parallel to the line of intersection of their planes. 

Since in the reflection from one mirror we are free to choose any normal as the 
axis of rotation, it is convenient to choose some point on the line of intersection of the 
two planes as the common foot of the normals to the two planes. With this choice, the 
equations prove the familiar result: The image space formed by successive reflections 
at two plane mirrors can be derived from the object space by rotating the object space 
through an angle 3 given by cos 40 = —cos (me, mi) about the line of intersection of the 
two planes. 

The effect of three mirrors is obtained by a further simple iteration: 
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i” = ngi’’ns = NgNoi’Nens = Ngh.Ninynon; = — gig’. (8) 


The third mirror introduces another central inversion. The negative sign before the 
last term in Eq. (8) represents the resultant central inversion. In this case we are free 
to choose the common point of intersection of the three planes as the common foot 
of the normals to the three planes. We may then state: The image space formed by suc- 
cessive reflections at three plane mirrors can be derived from the object space by rotating 
the object through an angle, 3 about an axis, a passing through the common point of inter- 
section of the three planes, followed by a central inversion about that point. The angle, 
0 and the direction of the axis, a of the resultant rotatory inversion are given by: 


cos $0 + asin 40 = g = ngnyn; = Snynn,; + Vn. (9) 
Now 
cos 40 = Sq = Sngynon; = — Sn; = Snonwn; = — Sngnwe 
= Snyngn, = — Smengn:. (10) 


If then, 40 is defined as a positive acute angle, one half the angle of rotation will 
have one of the four values, +48 and (180° + 48). The angle of rotation will therefore 
have one of the two values, +3. The magnitude of the angle of rotation is thus the 
same for all orders of incidence and its sign changes when the cyclic order is reversed. 

Synge’s convenient choice makes the subscript of each a correspond to the sub- 
script of the n of the second mirror on which the ray falls. Then by an elementary 
expansion formula: 


a; sin 408 = Vnonyn; = Vngnyn, = — n,Snem; + n2Sngh; + 3SniNo. (11) 
Letting 
M, = — Snon; = cos (Me, 03), M2 = — Snsn; = cos (Mz, 01) 
and M; = — Snim, = cos (m, m2), (12) 
we have 
— a, sin 38 = — Myn; + Mone + Mans, 
— a; sin 30 = + Mini — Mon, + Mans, (13) 


— a; sin 40 =-+ Mn; + Men: — M3ns3. 


These are obviously equivalent to Synge’s Eqs. (3.10) and show that in those equa- 
tions 


+ 1/k = sin 40. (14) 
Adding the second and third of Eqs. (13) gives 
— (a; + a3) sin 40 = 2M,n,, (15) 


so that the three vectors m;, a and a; are coplanar. Similarly, so are Ne, a; and a 
and ms, a: and a. 
The angles which n; makes with a and a; are given by 
sin 40 eos (m)@>) = — sin $9S5n)a. = M;, aaa M2M;3 —_- MoM; = M, 


= — sin $3Sn,a; = sin $0 cos (mas). 


(16) 
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o 
Thus a and a; are vectors coplanar with n; making the same angle with it. Since 
they are not identical they must lie on opposite sides of n. This proves that n, bisects 
the side aga; of the spherical triangle a,aa3, with similar results for nz and ns. 
Squaring the first of Eqs. (13) and remembering that a? =n? =n} =n? = —1, gives: 


— sin? 30 = — M; — Mz — Ms + 2M2M;Snon; — 2M;M,Snm; — 2MiM2Snymz, (17) 


Substituting from Eqs. (12) gives: 
sin? 48 = Mi + Mz + M3 — 2M:MoMs, (18) 


which is equivalent to Synge’s Eq. (3.9). 

Since, except for sign, ? is known from Eq. (18), the directions of the axes aj, a: 
and a; are, except for sense, given by Eqs. (13). The sense of each may be arbitrarily 
chosen. When it is chosen, the sign of sin $ for each axis is given. Combined with 
the sign of cos 40 determined from the order of incidence by Eqs. (10), this deter- 
mines the sign of # for each of the six orders. 

These formulae give in essence all that is contained in Synge’s three theorems, and 
the diagrams he gives can readily be constructed from them. Numerical computations 
using them are straightforward. They are especially simple when one or more of the 
angles between the mirror normals differ so little from 90° that small quantities of the 
second order may be neglected. 

In collimated systems, in which a group of any number of plane mirrors lies wholly 
within a region of parallel rays, any parallel displacement of the image space of the 
mirrors with respect to their object space does not affect the relative relation of the 
final image to the original object. For such systems the effect of any number of plane 
mirrors may be regarded as a rotation or rotatory inversion in which merely the angle 
of rotation and the direction of the axis are significant. These can also be found by 
simple iteration of Equation (4): 


if) = mypMy_1 ++ + MeMsiMyMe - - + Ny-M;, (19) 


giving 
(— 1)*¢ = (— 1)*(cos 438 + a sin $8) = Smumi-1- ++ MoM: + Vognyi-++ mem, (20) 


to determine the angle 8 and the direction of the axis of rotation, a of the rotation 
(Rk, even) or rotatory inversion (k, odd) of the object space into the image space after 
reflection from k mirrors. The computations are just as simple in form as those for 
three mirrors but of course the numerical work increases in proportion to the number 
of mirrors. 

An important special case of an even number of mirrors is a nondeviating erecting 
system, one which turns the image formed by a lens system through 180° without 
deviating rays parallel to the axis of the lens system. Neglecting, for the time being, 
the axial and lateral displacement of the image, the condition for this is that # = 180° 
or 

cos 30 = Sn.ne_1 ***HMon, = 0. (21) 


This may be written: ° 


Sn;(Snp_1 e+ Bish, a Vny_ oes N2N;) = S-n,Vmy-1 *** Hen, = 0. (22) 
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Since Vn,_1 - - - Mem, is parallel to the axis of the rotatory inversion produced by the 
first (k—1) mirrors, this merely requires that n; shall be perpendicular to that axis. 
We may therefore state: Any combination of an odd number of plane mirrors, with their 
normals arbitrarily oriented, may be converted to a nondeviating erecting system by adding 
another plane mirror whose normal is perpendicular to the axis of the rotatory inversion 
produced by the combination. 

An interesting, but impractical, special case is the combination of a cube corner 
with a fourth mirror in any orientation in which its normal makes an acute angle 
with the rays emerging from the cube corner. Since the quaternion product of three 
mutually perpendicular vectors is +1, it forms a nondeviating erecting system whose 
axis is parallel to the normal to the fourth mirror. 

The scalar part of the product of any number of quaternions (including vectors) 
is the same, so long as the cyclic order is maintained. Consequently the angle of rota- 
tion produced by any number of mirrors is the same for the same cyclic order of inci- 
dence, but the axis will in general be different for different orders. This permits stating 
a sufficient, but not necessary, special condition: 


O = Smyg-+-Mj°*+ My 
= Snin,Smz_-1 °°: > my* >> Be + S: Vnyn;,Vny_1 “0 Bae. > (23) 
a Sn jn j_15n j-2 + ° s+ ide * = * ees + S: Vnjn;_,Vn;_2 os ge > * * Bas 


This is obviously satisfied if each term on the right hand side of either of these equa- 
tions is equal to zero. We may therefore state: If, in a combination of an even number of 
plane mirrors, the normals to the first and last or to any two consecutive mirrors are per- 
pendicular to each other and the line of intersection of their planes is perpendicular to 
the resultant axis of the rotation which would be produced by the other mirrors alone, the 
combination will form.a nondeviating erecting system. 

The condition will obviously be satisfied if all the lines of intersection of the planes 
of the other mirrors are perpendicular to the line of intersection of the two. This even 
more special condition is met by all nondeviating erecting systems which I have found 
described. 

It is of course obvious that the mirrors must be so located as to permit unob- 
structed passage of the light through the combination. If a prism serves as a mirror 
combination it is further obvious that, to avoid refractive deviation, the exit face 
must be parallel to the image of the entrance face produced by the prism. This con- 
dition is of course never exactly met by any actual prism. One of the problems in prism 
design is to determine how large tolerances may be permitted in meeting such condi- 
tions without significant interference with the performance of an instrument. 

The effect of small deviations may be treated as the effect of thin glass wedges 
cemented to the surface of perfect prisms. Such a wedge will produce a rotation of the. 
image space on the object space represented by 

quXqu = (cos $3» + w sin 40,,)x(cos 0, — w sin #0,), (24) 
where the vector w is parallel to the edge of the wedge and #, is the angular deviation 
produced by it. Since these deviations are small it is usually sufficient to write 
Qu =(1+40,,w). This quaternion is inserted between the proper pair of normals in 
Eq. (19). The difference between the two results for 3, finite and 3, =0 gives the 
effect of the deviation. 
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Since, in general, four or more planes do not have a common point of intersection 
further computations are necessary if the relative position of image and object space 
is needed. 

To compute this we choose arbitrarily an origin of coordinates at d;_; =0 and de- 
fine object and image spaces by vectors r, (4=0, 1, 2,---,j,-:+,k—1, k) drawn 
from that origin. We represent by da a vector drawn from that origin to some arbi- 
trary point on the k-th mirro1. Equation (4) may then be written: 


i+ h = - thin ~ dhe . (25) 


Giving h in succession the values j, j7+1,---,k—1, k, (j=1,---,) and sub- 
stituting in succession the value of each ra_1 in Eq. (25) the effect of all or a part of 
the & mirrors is found. 

In carrying out these substitutions, it is convenient to use an abbreviated notation 
due to Hamilton. We define: 


ox = gxq-' = (cos 38 + a sin 30)x(cos 40 — a sin 30) (26) 


as an abbreviated form of Eqs. (1) and (2). Here ¢ is a “rotor,” a special form of linear 
vector operator which rotates a vector through an angle # about an axis parallel to 
the unit vector a. Carrying out the operations indicated in Eq. (26) and remembering 
that 

x = — aVax — aSax, (27) 
where —a Vax is the component of x perpendicular and —aSax, the component paral- 
lel to a, we find 


ox = (cos? + asin 3)(— aVax) + (— aSax) (28) 
x — ox = (1 — ¢)x = (1 — cosd — asin 3)(— aVax). ‘ (29) 


In the important special case when 3d =0, we have g=1, @x=x, and (1 —@)x obviously 
reduces to zero. Again, 


x + ox = (1+ ¢)x = (1+ cos? +a sin #)(— aVax) + 2(— aSax). (30) 


In the important special case when 3d = 180°, we have g= —1, and (1+¢@)x obviously 
reduces to 2(—aSax). We will later find it necessary to solve equations of the form 
of Eqs. (29) and (30) for the components of an unknown vector x. Since 


(1 — cos? + asin #)(1 — cos? — asin’) = 2(1 — cos¥), (31) 
Equation (29) is solved by: 


1 

aVax 7 ears (1 — cos? + asin d)(1 — ¢)x, (32) 
and the parallel component —aSax is obviously arbitrary. If } =0, this takes on the 
indeterminate form 0/0. However, as noted above, in that case, the operator, ¢=1, 
is a pure scalar and represents a parallel displacement as a rotation through a zero 
angle about an infinitely distant axis. Every vector in the object space is imaged by 
a parallel vector in the image space. Then, not only the parallel, but the perpendicular 
component of x is arbitrary. In this case the solution of Eq. (29) is wholly indeter- 
minate. 
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Equation (30) must first be solved for —aSax as follows. Remembering that 
SVax=0, SaVax=0 and a?= —1, we find 


— Sa(1 + ¢)x = — Sa(1 + cosd + asin 8)(— aVax) — 2Sa(— aSax) = — 2Sax. (33) 


Thus the component x parallel to a is one half, the component of (1+@)x parallel to 
the same axis. Then, analogous to the solution of Equation (29) we have 


1 
2(1 + cos #) 


If #=180° this takes on the indeterminate form 0/0. By returning to the original 
definition of ¢ in terms of gq, it is easy to show that in this case the perpendicular com- 
ponent —aVax is arbitrary and only the parallel component —aSax is defined by 
Eq. (30). The reflection from a single mirror represents a special case of d = 180°. In 
this case the vector a is normal to the mirror: Adding to x any vector in the plane of 
the mirror merely displaces the foot of the normal which, as noted at the beginning, is 
arbitrary. If ? 180°, both components are determinate and their sum gives the vec- 
tor x. It is readily seen that: 


(1 + cos # — asin #)[(1 + ¢)x + aSa(1+¢)x]. (34) 


— aVax = 


o(x + y) = ox + oy (35) 
and 
-1, —1 —1 -—1 . 
PiPrX = Qi(GaXGn Qi = GiGrXQn Qi - (36) 
The product of several such operations in sequence is thus evaluated by the multi- 


plication in proper order of the quaternions in terms of which they are defined. We 
may then define: 


on -—1 —1 -1 -l 
or iX = qniXQn + = MpMj1* + * Mig:MsXN; Nigy- ++ MaawM, . (37) 
It is frequently convenient to note: ; 
Ph, iP ci-1),5 = Oh, i- (38) 
With this notation Eq. (25) may be written: 
Tn = — Gr wEr—1 + (1 + ona) da. (39) 
Starting with 
tj = — ;,2j3-1 + (1 + 6;,;)d;, (40) 
and substituting successively the values of rj, fj41, - + + , Te-1 for fr in Eq. (39), we 
find 
hak 
re = (— 1) gy ia + DO [1 + (— 1)* bea] (da — da). (41) 
h=j 


If 7=1 and all the d,’s are zero this obviously reduces to Eq. (19). 
The process of computation consists in calculating in order 


Qkk = Mn,°°* » Ven = Qk r4iMan, *** » Qk, = Qk,2Mh1, (42) 


from which each of the 3;,,’s and a%,,’s are computed by: 
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Sqn = COS $0x,n, arn = Vaen/sin 30x,n. (43) 


The computations may frequently be shortened by using Eq. (38). The term @«,j-1 
is then computed either by Eq. (26) or by Eq. (28) whichever seems easier, and each 
of the terms in the summation in Eq. (41) by Eq. (29) or (30). They are then added. 

Since the origin of coordinates and the points to which the vectors d, are drawn 
are arbitrary, the computations can ordinarily be much simplified by suitably choos- 
ing them. 

When two successive mirrors intersect, choosing a point on their common inter- 
section makes d,—d,-1=0 and the corresponding term in the summation vanishes. 
When three successive mirrors intersect, choosing their common point of intersection 
as the common foot of their normals makes d,_1 = d,=da41 and two successive terms 
of the summation vanish. If they are the first three and we choose it as the origin of 
coordinates, we have d;.1.=d;=dj,:=dj,2.=0 and three successive terms vanish. It 
was this choice which made the treatment of the three mirror problem so much 
simpler. 

In numerical work it is convenient to express all the vectors in terms of a suitably 
chosen rectangular coordinate system represented by unit vectors, i, j, and k. By 
substituting for rj, i, j, and k in turn and carrying through the computation, the 
position of the image of any point, rj.1=xi+yj+<zk, in the object space may be ob- 
tained by simple addition. 

Equations (28), (29), (30), (38), (41), (42), and (43) are thus all that are needed 
to determine the relative position and orientation of object and image after reflection 
from any number of plane mirrors. It is, however, sometimes more convenient to 
calculate the parameters of the transformation of the object space into the image 
space. 

The most general rigid displacement of a three dimensional Euclidean configura- 
tion is a twist about a screw. The relative position of image and object after reflection 
from an even number of plane mirrors will be known when the angle of the twist and 
the pitch and the position of the axis of the resultant screw are known. 

A twist through an angle 3; about a screw of pitch a,/9; whose axis is parallel to 
the unit vector a; is equivalent to a rotation about that axis followed by a parallel 
displacement a,a; along the axis. The effect of an even number of mirrors will there- 
fore be given by an equation of the form: 


Te = PiFj-1 + (1 — oe) de, + ora, (44) 
where the vector, d;,; from.the origin to the axis of the twist is given the double sub- 
script to distinguish it from the vector d, drawn to a point on the k-th mirror. Com- 
parison with Eq. (41) shows that 


de = Ok,j, Be = Vx,; and a= Qk, ;- (45) 
These are computed by Eqs. (42) and (43) and a, is given by 
h=k 
a= — Sa. >, f1+(- 1)*-"og.n] (da — d,_1). (46) 
h=j 
Then 
h=k 
(1 — px)dk,; [1 + (— 1)**pan](da — das) — nae. (47) 


h=j 
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These equations are solved for the perpendicular components of the dx,;’s, by Eq. (32). 
Their parallel components are arbitrary. The effect of all the even combinations of 
the k mirrors is then known in terms of the angle of the twist and the pitch and the 
position of the axis of the resultant screws. 

If a twist upon a screw is combined with a central inversion, the result is a rotatory 
inversion. The relative position of image and object after reflection from an odd num- 
ber of mirrors will be known when the position of the axis, the location of the center 
on that axis, and the angle of rotation of the resultant rotatory inversion are known. 

A rotatory inversion through an angle 3; about an axis parallel to a unit vector a; 
is readily shown to be given by an equation of the form 


Te = — Oitj-1 + (1 + ed dk,;, (48) 


where d,,; is a vector drawn from the origin to the center of the resultant inversion. 
Comparison with Eq. (41) shows that, as in the cast of an even number of mirrors. 


od, = Pk, is 0, = Ok, j and ak = ax,j- (49) 
These are computed by Eqs. (42) and (43). Then 


h=k 
(1 + pi)de,j = Do [1 + (— 1)* hen] (da — das). (50) 
h=j 


These equations are solved for the parallel and perpendicular components of the 
d,,;'s by Eqs. (33) and (34). The effect of all the odd combinations of the k mirrors 








Fic. 1. Range finder five mirror prism combination. 


is then known in terms of the position of the axis, the location of the center of inver- 
sion on that axis and the angle of rotation of their resultant rotatory inversions. 

Although there are necessarily a considerable number of steps in such a computa- 
tion, the steps are not complicated, consisting merely of simple multiplications and 
separation of quaternions into their scalar and vector parts. 

As an illustration of the straightforwardness of these quaternion computations, 
suitable dimensions for the five-mirror combination shown in Fig. 1 are computed 
below. This is part of the prism system used in an Army range finder. The other part 
is a simple two-mirror system cemented below Face 1. 

In use the axial ray from the right hand objective enters through Face 6, is re- 
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flected in turn from Faces 5, 4, 3 and 2 in succession and from the middle of the near 
edge of the silver strip on Face 1 at the point F in Fig. 1, and finally emerges through 
Face 2, which thus serves both as a mirror and an exit face. The image is focussed on 
that edge of the silver strip, where it is brought into coincidence with the image from 
the left hand objective, coming through the lower two-mirror system which is not 
shown. 

We must assure ourselves that the image of the exit Face 2, in the space of the 
objective, is parallel to the entrance Face 6 and that light from the whole aperture 
of the objective has unobstructed passage to the edge of the silver strip where the 
coincidence is observed, and for a reasonable distance around it. 

We therefore compute, in the space of the objective, the direction of the image 
of the normal to Face 2 and the position of the corners of the silver strip and of any 
prism faces which might obstruct the light. The plane of a mirror and its image in 
that plane coincide, so that we need only compute the relationship between the space 
of the objective and every other object space in which the prism faces lie. 

Since the planes of Faces 1, 2 and 3 produced, intersect in a point, it is convenient 
to choose that point as the origin of co-ordinates. Choosing axes i, j, and k and di- 
mensions as shown, and remembering that cos 60°=sin 30°=}3, sin 60°=cos 30° 


=4/3/2, and cos 45°=sin 45°=1//2, we have: 


as 1 
ni = k, d; = 0; No = $(j mm 1 3k), ds — 0; n3; = /9 (i = j); d; = (0; 
Vv 


i 1 
n,=- get e, d, = (a + 2b)i + hk; ns = +7 (i+k), ds = (a + 26)i, 
V V 


where h=a tan 30° =a/v/3. 
Then from Eq. (41), 
Ts — $5,403 = 5. + 5,882 = T5 + 5,:To = (1 + 5,5)(ds — da) + (1 — 5,4)(da — ds). 


By Eqs. (42) and (43), 


1 1 
9s5=N5=—=(i+ k), cos }d55=0, 35,5 = 180°, a55=n5 = —(it+ k). 
v2 v2 


Since #;,,= 180°, we have by Eq. (30). 


1 1 
(1 + o5,5) (ds - d,) 2[- @5,5585,5(d5 — d,) | = | - oman (i + k)S : (i +- k(t | 
VJ2 V2 


— hi+k). 


It may be interesting to note again here that, since the position of the foot of the nor- 
mal is arbitrary, d, and ds could have been written more generally as d4=(a+26 
+h—z)i+yj+zk, ds=(a+2b—2z’)i+y’j+z2’k, and (1+¢5,5)(ds—d,) would, as the re- 
sult of more complicated computations, have come out the same. 


95.4 = 95,04 = — 3(i+k)\(i+k) = 1. 





a ee 





f 


— 
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Then according to Eq. (29), ds4fs=fs and (1 —@5,4) (ds —ds3) =O and Eq. (41) reduces 
to 


T5 = 3 — h(i + k) = — $5,3%2 — Ali + k) = — $5,1To — h(i + k). 


1 
%3 = a3 = — (i — j), 


J/2 
By Eq. (26), 
— $5,3(X2i + yoj + 22k) 


B(i + j)(xei + yoj + 22k) (i — j) = yoi + xej + ook 
Pee ae ree EP ae 
1 = ganas = 7 — HG vse = 1+ V3i V3j+k) 
— 4-14 V3i — V3 j + k)(xoi + yoj + zok) 
(-1-V3i+V3j-&) 
4(yo a V3 Z0)i + Xoj + 3(z0 + V3 yo)k. 


A unit vector normal to the exit Face 2 of the prism is n:=4}(j—+/3k) giving x»=0, 


— 5,1(Xoi + voj + zok) 


yo=} and 29= —+/3/2. Its image in the space of the objective is therefore 
ifi _({ v3 i a 

— $5,122 = —| — — Vz — — i+oj+—[-~= —_ 3 ie = i, 
— TF vi >) wrt ee 


which is normal to the entrance Face 6 of the prism. If the angles of the prism com- 
bination are correct there will thus be no prismatic deviation. 
The corners of the silver strip also lie in the space of ro. Their images are there- 


fore given by 
— $s,:[ai + (c+ w)j] — htit+), — $5,:[(@ + 2b)i + (c + w)j] — h(i+k), 
— $5,:[ai + cj] —hit+h), — $5,1[(a + 2b)i + cj] — hi+k). 
Carrying out the computations we find 
Ac + w — 2h)i + aj + $[V3 (c + w) — 2h)k, 
3(c + w — 2h)i + (a + 26)j + [V3 (c + w) — 2h)k, 
. 3(c — 2h)i + aj + [V3 ¢ — 2h]k, 
Ac — 2h)it+ (a + 26)j + 4[V3c — 2h]k. 
These give the location, in the space of the objective, of the “air image” of the silver 


strip, i.e., the location the image would have if the prism faces were front surface 


mirrors in air. 
Face 5 lies in the space of the objective. Its corners and their “air images” are 


given by 
(a + 2b)i + aj, (a + 26)i + (a + 2)j, 
(a + 2b + e)i + aj — ek, (a + 2b + e)i+ (a + 2b)j — ek. 


If the diameter d, of the objective is greater than 2b, which is the case in all prac- 
tical constructions, and light from the whole aperture of the objective is to reach all 
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parts of the silver strip, it is obvious, since the light is converging to a focus, that the 
j and k co-ordinates of the image of the silver strip must lie somewhat within the 
limits of the corresponding corners of Face 5. Since their j co-ordinates are equal it 
is obvious that there will necessarily be some vignetting of the image at the ends of 
the silver strip, the amount depending upon the aperture and focal length of the 
objective and the size of the prisms. 

For the k co-ordinate of the back edge of this strip we have: 


41/3 (c + w) — 2h] < 0 


or 
(¢ + w) < 2hk/V3 = 20/3, 

Allowing for some vignetting at the rear edge of the silver strip, which will do no 
harm, we may set (c+w) =2a/3. To allow an equal area for the matching image from 
the other objective we may set c=a/3. 

For economy of glass it is obviously desirable to make the viewing aperture ap- 
proximately square. The projection of (c+w) =2a/3 on Face 2 is »/3-2a, 3=a/V3. 
For economy of glass we should therefore make b=a/2+/3. 

The bottom edge of Face 5 coincides with the bottom edge of Face 6. The image is 
focussed on the bottom edge of the silver strip. The length of the path of the axial 
ray.in the glass is therefore the difference of their i co-ordinates: 


(a + 2b +) — 4c — 2h) = (5 + 4V3)a/6 4+ €. 


Because of the index of refraction yu of the glass the image of the focal plane as 
seen from the lens in air will be displaced in the +idirection by (u—1)/u times this 
length. The effective image of the focal plane, in the space of the objective, lies there- 
fore at a distance 


[1 — (u — 1)/u][(5 + 4V/3)a/6 + e] = [(5 + 40/3) 0/6 + e]/u 


in the —i direction back of Face 6. The back principal plane of the objective must be 
placed a distance equal to its focal length f further in the +i direction. 

At the front edge of the silver strip, where the images coincide, vignetting would 
be intolerable, so we must have 


4[V/3c— 2h] >—-—e . 


or 
e>h—vV3c/2 = h/2 = a/2V3 = V3 0/6, 

The thickness e of the rhomb forming Faces 4 and 5 of the mirror combination must, 

therefore, be somewhat greater than one half the height of the rest of the combination. 
The convergence of the rays in the glass is determined by the distance between 

the focal plane and the “glass image” of the lens, i.e., the distance the back principal 

plane would have if the light path to it lay wholly in glass. This distance is yf. If d 

is the diameter of the lens the theoretically necessary value is given by; 


e = V3a/6 + [(5 + 4V/3)a/6 + e]d/2uf 


or 
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e = [V3 + (5 + 4V/3)d/2uf]a/6(1 — d/2y/f). 


In practice e is made slightly larger to provide against small misalignments. 

In this prism combination, computation of the location of the images of the 
corners of the other prism faces is not necessary, because simple geometrical consid- 
erations show that Face 5 forms the limiting aperture. However, they could readily 
be computed if desired. 

To complete the design it would, of course, be necessary to compute the dimen- 
sions and shape of the lower two-mirror prism so that the axial ray from the left hand 
objective lies in the same line, 


xi + (a + b)j + (V3 c/2 — A)k, 


as the axial ray from the right hand objective, after two reflections in the prism has 






































Fic. 2. Leman-Sprenger nondeviating 
erecting prism. 


the direction, — nm, and is focussed on the same point F, giving r? = (a+b)i+cj. These 
computations would obviously be much shorter than the previous ones. 

For the purpose of establishing tolerances, computation of the effect of small 
deviations of the faces from their correct orientation is necessarily longer and more 
tedious, but is just as straightforward. 

These computations of suitable dimensions of this prism combination could, of 
course, have been made by determining the parameters, the axis and angle of the 
rotatory inversion produced by Faces 1, 2 and 3, using Eqs. (10), (13) and (18), and 
combining their effect with the parallel displacement produced by Faces 4 and 5. The 
computations would, however, have been much longer. 
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As an illustration of the computation of the parameters of a transformation of an 
object space into an image space, the effect of a Leman-Sprenger nondeviating erect- 
ing prism is calculated below. Axes i, j and k and distances are chosen as shown in 
Fig. 2. Remembering that: 


cos 60° = sin 30° = 1/2, sin 60° = cos 30° = 1/3/2 and 


sin 45° = cos 45° = 1/)/2, 


we have: 
fm @ 1 . . . 
n, = }(V3i—k), d, = 0; as = 5 (—-it dD), d; = 4i — j, 
Vv 
a: 1 4 . . . 
nm = 3(V3it+k), de = 0; Bee lm io Di d, = 4i — j. 
V 


The combination of Eqs. (41) and (44) reduces to: 
I, — Glo = 14 — G4,:To = (1 — Ga,3)(ds — de) = (1 — da)dayr + oats. 


By Eq. (37), 


Qus = OM; = 3(—i— j)(—it+j) = —k, 
92,1 = Moms = 2(/3i1+ k)(V/3i — k) = 3(-—14+ V3 5). 
By Eq. (38) 
gs = 94,1 = 94,2921 = — $k(— 1+ V3j) = #(V3i+ bk), 
giving: 
cos 330, = — Sq, = 0, cos $043 = — Sqas = 0, 


1, cos }4 = cosd,3 = — l, sin 34 = sin 34,3 = 0. 


sin 434 = sin 304.3 


By Eq. (43), 


a,= Vq,/sin 30, = 3(y 3i + k), ag3 = Vqa,3/sin 4043 = — k, 
Then ’ 
— a4,3Va4,3(ds — de) = — (— k)V(— k)(4i — j) = k(— 4j — i) = (4i — jp. 


By Eq. (29), 
(1 — $4,s)(ds — de) = (1 — cos d4,3 + a4,3 Sin 34,3)(4i — j) 
2(4i — j) = 8i — 2j = (1 — da) dar + oats. 


By Eq. (46), 
a, = — Sag(1 — $4,s)(ds — de) = $S(\/3 i + k)(8i — 2j) = 4V3. 


By Eq. (47), 


(1 — a)da,1 = (1 — Ga,3)(ds — de) — aya, 
= (81 — 2j) — 2/3 (V3i+k) = 2i- j — V3 k). 
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By Eq. (32), 


~ (1 — Cos Da + a, sin 8,)(1 -~ 4) 41 =i-— j = 2/3 k. 


—_ a,Vasda 1 SE Sean: 
2(1 — cos &4) 


Since the component of d4, parallel to the axis a, is arbitrary it is convenient to 
add to this normal component a parallel component 


2/3 a, = 3i+ 2V3 k, 
giving for dy, the value 
di = 4i —.j. 


The resultant axis of rotation, shown in Fig. 2, thus intersects the prolongation of the 
roof angle formed by ng and nm, at the point r=4i—j and is perpendicular to the en- 


tance face of the prism. 
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Fic. 3. Leman-Sprenger nondeviating erecting prism. 


If the erecting system consisted of front surface mirrors, the image, besides being 
rotated through 180° about the axis would be displaced backward through the dis- 
tance a,=4,/3. Since, however, the path through the prism is in glass, this backward 
displacement is partially compensated for by a forward displacement equal to 
(u—1)/p times the length of the glass path, where yp is the index of refraction of the 
glass. The glass path is readily determined by rotating the entrance face AB (Fig. 3) 
about the axis to A’B’ and displacing it backward through the distance a= 4/3 
to A’’B’’. The distance of this “air image” of the entrance face from the exit 
face CD is the length of the glass path and is readily seen to be 6/3. Assuming 


w=1.5, 6V/3(u— 1)/u= 24/3. The relative positions of object, “air image” and image 


are shown in Fig. 3. 
The effect of a perfect prism of this type could obviously be determined much 
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more easily by a simple geometrical construction. Further, the quaternion computa- 
tions have obviously been carried through in needless detail for a perfect prism. This 
has been done to show the steps which would be necessary for the purpose of estab- 
lishing tolerances in determining the effect of small deviations of the faces from their 
correct orientation. In such computations, the scalar part of the quaternions g4,3 and 
gs would not be zero, but small quantities representing the effect of the deviation of 
the prism faces from their proper orientation, making necessary the computation of 
the perpendicular and parallel components of vectors which, in the example given, 
could have been written down by inspection. 
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ROTOR BLADE FLAPPING MOTION* 


BY 


GABRIEL HORVAY 
McDonnell Aircraft Corporation 


1. Introduction. In an earlier paper a method was developed for the solution of a 
class of Hill differential equations.' In the present paper the method is applied to the 
solution of the differential equation of flapping motion of helicopter rotor blades. 
A brief review of the equation in Sec. 2 is followed by the solution of the homogeneous 
equation in Sec. 3. In Sec. 4 the stability of blade flapping motion is established for 
the range of the numerical parameters which usually arise in helicopter theory. The 
stability of the equation for large values of yu is discussed in Sec. 5. In Sec. 6 the possi- 
bility of resonance is investigated. 

2. Differential equation of rotor blade flapping motion. The blades of a helicopter 
rotor are not attached rigidly to the rotor shaft, but are hinged so as to permit the 
blades to flap out of the rotor plane. This prevents transmission of severe alternating 
moments from the blades to the rotor shaft. The angle between a blade and the rotor 
plane (the plane normal to the rotor shaft) is called the flapping angle 6; it is counted 
positive upwards. 

At a standstill, the blades of the rotor rest on the droop stops. When the rotor is 
revolving at the angular velocity w, the centrifugal force 


R R 
f dC = f rw*dm (1a) 
0 0 


(dm is the mass of the blade element at distance r from the axis of rotation, R the 
blade length) lifts the blades off the stops and keeps them in approximately horizontal 
position. The blades also experience lift forces 


R 1 fF 
f aL = ~f pca(y + 8) U 7dr. (1b) 
0 0 


Here p is the air density, c the chord of the blade profile, a the slope of the lift curve, 
3% the blade pitch angle, g¢ + the angle of attack, and U the resultant air velocity, 
cf. Fig. 1a. 

In vertical flight the equilibrium angle 8 of the blade is determined from the condi- 
tion, cf. Fig. 1b, that the moments of lift, weight and centrifugal forces balance each 
other: 


* Received Sept. 5, 1946. This paper is based on an investigation which was carried out at the 
McDonnell Aircraft Corporation under sponsorship of the Bureau of Aeronautics, U. S. Navy Depart- 
ment. 

1G. Horvay, Unstable solutions of a class of Hill differential equations, Quarterly Appl. Math. 000, 
385-396 (1947). In the following, this paper will be referred to as I. (Note that A in Eg. (29) of I should 
be written with a superscript 2 instead of 0.) The present paper isa continuation of the study, and deals 
with the practical aspects of the problem. The writer's thanks are due to his colleagues, Elizabeth J. 
Spitzer, Kathryn Meyer, and Frances Schmitz for checking the derivations and for help with the numeri- 


cal calculations. 
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R 
f r(dL — gdm — BdC) = 0. 
0 


(The drag has a negligible effect on the flapping motion.) It is assumed here as else- 
where that the blade is rigid, and that 8 is so small that the approximation sin B~8, 
cos 81 is valid. The same assumption is made also for the angles # and ¢.? 

Let ¢ denote the time and 


¥ = ot (2) 


the azimuth of the blade as measured from the rearmost position. In the course of 
forward flight at the speed of wwR (yu, the so-called advance ratio, is the ratio of forward 


AL 











U 





(a) (¢) 


Fic. 1. Velocities and forces at a blade element. 


velocity to blade tip speed), an advancing blade element (¥=90°) moves with the 
speed w (r+uR), a retreating element (~=270°) moves with the speed w(r—yR). The 
variation in the relative air velocity brings forth a variable lift moment, and thus 
produces a blade flapping which varies with the azimuth y. The tangential com- 
ponent of the air velocity at element dm is 


Ur = rw + poR sin y, (3a) 


and the perpendicular component is 
Up = wR — 18 — woRB cosy. (3b) 
Here \wR is the difference between the sinking speed of the helicopter and the induced 


velocity through the rotor disc. \ is called the inflow ratio. 
The flapping motion of the blade is governed by the equation 


R 
Ip = f r(dL — gdm — pdC), (4) 
0 


where 


R 
I= f r?dm (4a) 
0 


is the moment of inertia of the blade with respect to the flapping hinge. Observing that 


2 A brief list of the assumptions involved is given at the end of Section 4. For a more extended dis- 
cussion cf. J. B. Wheatley, Aerodynamic analysis of the autogiro rotor, NACA Technical Report 487, 1934. 
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(o + 8)U2~ UrUp + 8Ur (3c) 


performing the integration (4a), and dividing by J, one arrives at the differential 
equation of flapping motion 


F(8) = E(t), (5) 
where 
F(8) =B + pb + s(B (Sa) 
with 
p(t) = no[1+4usin wt], s(t) = w?[1 + $n cos wt + mu? sin 2wt]; (Sb) 
and : 
E(t) = — mgre,/I + nw®{4r + (1 + w2)0 + (2Au + Sud) sin wt — wd cos 2wt}. (5c) 


The first term of the forcing function E(wt) represents the moment of the blade weight 
(r-g is the c.g. distance from the hinge), and is usually negligibly small. The remaining 
terms constitute the aerodynamic excitation. Control of the helicopter is effected by 
cyclic pitch variation: 

3 = do + BJ, cos wi + IV, sin wt. (5d) 


Here @ is the “collective pitch” setting, 3. the lateral, —d, the fore-and-aft control 
setting. 
It will be convenient to call the equation 


F(8) = 0 (6) 


the “homogenized” flapping equation. The term B of F (8) is the inertia term. The ex- 
pression p(t)8 represents the damping which is produced by the change insthe angle of 
attack caused by flapping. The dimensionless quantity 


n = pcaR*/8I (7) 


can be called the aerodynamic damping coefficient. The leading member, wf, of the 
spring term s(t)8 represents the restoring moment of the centrifugal force. The fre- 
quency w of the coefficients p(t) and s(t) can be called the parametric frequency.‘ 
In the present problem this frequency coincides with the frequency impressed by the 
function E(w). 

The solution of Eq. (5) can be written in the form 


B = Bo + AB: + heBe (8) 


where {> is the particular integral and represents the forced (steady state) response; 
8, and @ are the two solutions of the homogeneous equation, and represent the natu- 
ral modes (transients). 4; and he are numerical coefficients determined by the initial 
conditions. 

Two papers appear in the literature on the solution of the homogeneous equation, 


® The Germans refer to y=8n as the “Locksche Trigheitszahl.” 
‘ N. Minorsky, On parametric excitation, Journal Franklin Institute 250, 25, 1945. 








/ 
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7(8) =0. Bennett® has investigated the solution for the particular values n=1.5, 
u=1.0, B(0) =1, B(0) =0, in the following manner. He determined the coefficients of 


the expansion 
B(Y) = B(0) + 48(0)/1! + #8(0)/21+ - -- (9) 


by repeated differentiation of Eq. (6), plotted B(W) from y=0 to y=4r, and noted 
that 8(47)/8(0) <1. From this result Bennett concluded that stability is insured for 
all 1<1. In an earlier paper Glauert and Shone® solved Eq. (6) by omitting the 8 
term for simplicity, and concluded that the motion is unstable. 

Evidently, the above approaches to the stability problem are approximate and 
inconclusive. 

More extended investigations are available for the particular integral. By can be 
assumed in the form’ 


Bo = do + a; cos wt + Bb; sin wt + a2 cos 2wt + be sin 2wi+---. (10) 


The constant part, do, is called the coning angle of the blade. The angles a; and 3}, 
determine the tilt of the cone axis (forward and to the port side, respectively, when w 
is counterclockwise) ; the higher terms determine the motion of the blade in and out 
of the cone. Substituting (10) into Eq. (5), one obtains the infinite system of equations 


1: ado + dp?nbe = — megr.g/Tw? + $ndrd + (1 + wu?) ndo + gund, 
cosW: $udo + (1 + 3u?)bi — Fuae + 3u7bs = (1 + fu?)d. 
siny: — (1 — $u*)a1 — Jub, — fu7a3 = 2wrd + Sudo + (1 + Fu’)d, 
cos 2y: una, — 3a, + 2nbe — $unaz3 +--- = — w?ndo — $und, (11) 
sin 2: p?ndo + $unb; — 2naz — 3b2 — $unbs + --- = fund, 
cos 3p: 4u2nb; — 2una, + 8a3 — 3nbs+--- = — Zu?nd, 
sin 3y: 4yu?na; + 2unbe — 3na3 — 8b3 + --: = — Su?nd, 
in the flapping coefficients do, a1, bi, - - - . Placing a, =b,=0 for n>WN, one can solve 
the first 2NV+1 equations for do, a:, - - « , av, by. For instance, for the numerical values 
n = 1.7, megre,/Iw? = 0.03, } = — 0.10, do = 0.2, 8. = 9, = 0, w = 0.34738 (12a) 


which are representative for a helicopter, one finds that 
Bo(W) = 0.124 — 0.125 cos Y — 0.057 sin y — 0.012 cos 2y 
+ 0.007 sin 2 — 0.001 cos 3y +: -- (12b) 


is the forced response. 
Glauert® initiated the above method of approach in his fundamental paper on 


5 J. A.J. Bennett, Rotary-wing aircraft, Aircraft Engineering, May, 1940. An account of the work of 
Glauert and Shone‘ is also given in this paper. 

6 GJauert and Shone, R.A.E. Report BA 1080, 1933. 

7 This is a departure from the commonly accepted notation Bp =a@o—a; Cos wt —b; sin wt—az cos 2wt 
—besin2wt—---. 

® H. A. Glauert, A general theory of the autogyro, R. & M. No. 1111, British A.R.C., 1926. 
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rotor blade theory, and carried the expansion of Bo to first harmonics.- Subsequently 
Lock,® Wheatly? and others determined also the second and higher harmonic flapping 
coefficients. They found from numerous examples (for 450.5) that the convergence 
of 8» is sufficiently rapid to permit termination of the series with the second harmonic 
terms. 

Although a numerical determination of By) has thus been accomplished, (in con- 
trast with the lack of similar solutions for 8; and $2,) the problem of forced vibrations 
cannot be considered as completely solved so long as the convergence of series (10) 
is based solely on an appeal to numerical examples. It will be shown in Section 6 that 
for definite values of » and m the expansion (10) diverges. The motion B is then said 
to be in resonance with the excitation E(t). 

3. Solution of the homogeneous equation. Determination of the transients is very 
simple in vertical flight. Placing »=0 in Eq. (6) one obtains the familiar equation of 


damped motion 


d*p/dy? + ndB/dy + B = 0. (13) 
This has the solutions 
Bi = exp (— idm vi- i0¥, (14) 
B, = exp (— 3np)-sin V1 — jn? y 


For 4 >0 one can write in analogy 
= P : 
Bi e:(y) i(y) (15a, b) 
Bo = e2(p)P2(y), 
where 
ex(y) = exp (— 3Napp¥); (16a, b) 
e2(y) exp (— 3Mapps¥); 


and P;(W) and P2() are oscillatory functions. map, can be called the apparent damping 


coefficient. 
The functions @, é2, P:, P2 can be determined in the following manner. One sub- 


10,11 


stitutes 


By) = o(f) -e(y), (17) 
1 t 
e(y) = exp (- —f p(bdt) = exp (— nv + §nu cos y) (17a) 


into Eq. (6). This leads to an equation for the factor v(y) in which the first derivative 
is absent. Replacing the trigonometric functions by exponentials, the equation can 
be written in the form 


d?v/dy? + [dre + Oe + Oy + Ore + Boe*¥]y = 0, (18) 


°C. N.H. Lock, Further development of autogyro theory, R. & M. No. 1127, British A.R.C. 1928. 
10 Whittaker and Watson, Modern analysis, p. 194, Cambridge, 1927. 
1 The notation {* indicates that only the upper limit, ¢, is to be substituted into the indefinite in- 


tegral of p(t). 
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where 
0 = 1 — jn? — gy’n’, 
6, = Zun(1 + in), (19a, b, c) 
6. = u*n(4n — 41), 
and 6* denotes the conjugate complex of 9. By writing v(W) in the form 
+00 
v(y) = e% > Cret*¥ (20) 


the differential equation (18) is reduced to an infinite system of linear homogeneous 
equations in the C,. 
For the values 
i q ; 
5 o = — arctan >So + mi for —1 =< q 
T / 1 — @’ 


lA 


(21a) 


1 en 
—log (q+ Vg? — 1) + (m—})i for gS—1,¢g21 (21b) 
T 


= * log (q + Vq’2? +1) + mi for q imaginary (21c) 
where 
m=0,+1,+2,::- (22a) 
g=VDsinrV/6 or VW—Dsinh rV/— 6 (22b) 
g =V—DsinrV% or VWDsinh r/— (22c) 
ye = 1/(00 — k?) (22d) 
and!” 
1 Oi*y2 Oo*y2 =O 0 
A191 1 O:*y1 O*y, O 
D=|- B2y0 91¥0 1 O;*yo Ostyo = (|, (22e) 
0 Oy: Oy 1 Oty, 
0 0 B22 Ai ¥2 1 








the system of equations is consistent, and the expansion coefficients C;, can be de- 
termined. Two solutions are obtained, Cf and C. The first set for, say, 


o = 6; = o, + 190i, o- 2 0, (23a) 


2 It was shown in I that the infinite determinant DD may be expanded (for small 4) into a power 
series in 6=| 6,| 2, e= | 62|2 and 41=}(0102+6; 02). The leading term of the expansion is 1; the higher 
coefficients are functions of @) only and may be obtained from the Tables of I. 











1947} ROTOR BLADE FLAPPING MOTION 155 


the second set for 
o¢ = o2= — 61. (23b) 


The two solutions are linearly independent, except when ¢,=0, and 20;=integer. 
The trigonometric function P;(W) is now obtained as the product 


+s 
Pi(y) = exp (3ny cos y)- : ce” exp i(o; + k)y (24a) 


k=—oo 
and the exponential function e,(W) as the product 
ex(¥) = exp (— gmp) -exp (¢-¥). (24b) 
The determination of Pe(W) and e2(W) (call these Eqs. (25a, b)) is similar. It is seen 
now that e; involves the apparent damping coefficient 


Rapp = 0% — 20, (26a) 


which is smaller than n, and é involves the apparent damping coefficient 
Napp, = n + 20, (26b) 


which is larger than n. Evidently, map», is the critical damping coefficient.” It will be 
written without the subscript 1. 

Example. Determine the flapping transients 6; and 62 of a helicopter for which 
n=1.7. Choose an advance ratio u between 0.30 and 0.35 so that the Tables of I 
(which cover the range 6) = —0.6 to +0.5 in 0.05 steps) be directly applicable. 

One finds from (19a) that for 


n = 1.7 (27a) 
and 
0 = 0.2 (27b) 
the advance ratio has the value 
uw = 0.34738. (27c) 
This yields, by (19b, c) the values 
6, = 0.19685 + 0.334657, 6. = 0.03875 — 0.10258:. (27d) 
The characteristic exponents 
o, = 0.30782 + 41, = — 0 (28) 


associated with the numerical values (27), were determined in I, Section 4. The cor- 
responding functions v,(W) and v2(W) were also given there. Eqs. (24) and (25) now 
yield 





8 » —2¢,<0 implies instability. The arbitrariness in the value of m, Eq. (22a), is removed by 
making Co the dominant term of v(y). 
“4 The normalization is used: the sin 4y term of v; and the cos 4y term of vz are written with the coeffi- 


cient 1. 
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€1(y) = e~0-5422y, (29a) 
Pi(y) = — 0.0733 cos 4y + 0.8369 sin 4y + 0.1923 cos fy + 0.1717 sin $y 
+ 0.0268 cos gy + 0.0131 sin $y + 0.0023 cos $y + 0.0023 sin gy 
+ 0.0002 cos gy + 0.0004 sin $y +---, (29b) 
and 
eo(y) = e7!-15789, (29c) 
P(y) = 1.2863 cos 4y + 0.3795 sin $y + 0.4404 cos $y + 0.1028 sin $y 


+ 0.0676 cos §¥ + 0.0229 sin $y + 0.0068 cos $y + 0.0050 sin fy 
+ 0.0005 cos $y + 0.0006 sin$y+---. (29d) 
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Fic. 2. Transient blade flapping modes 6:(y) =e:(y)Pily) and B2(y) =e2(¥)P2(y) for 
damping coefficient nm =1.7 and advance ratio u =0.34738. 
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The factors ¢:, P; and é:, Pz of the transient flapping angles 6; and {2 are plotted 
in Fig. 2 as functions of y. Pi(y) and P(y) are drawn normalized to 1 at their maxima, 
e:(W) and e:(W) are drawn normalized to 1 at ¥y = 360°. The graphs indicate that at the 
very large advance ratio 4 =0.347 the flapping motion of the blade, as given by Eq. 
(6), is still so stable that a transient reduces to less than 4% of its original amplitude 
within one revolution. The wave form P(W) of the transient contains a pronounced 
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third harmonic component and unimportant higher harmonics. The frequency of the 
transient (or more exactly, of the dominant term of the transient) is one half of the 
rotor speed. 

4. Stability of blade flapping motion. It was seen in the preceding section that the 
transient f; is associated with an apparent damping coefficient mapy, which is less than 
n when o, 0. o, can be called the “absolute degree of destabilization” of the system. 
For a damped system (n>0), like the present one, it is more convenient to charac- 
terize destabilization in terms of the ratio 


2o,/n = 1 — mapp/N. (30) 


This can be called the “relative degree of destabilization” of the system, or briefly, 
the degree of destabilization. For the example of the preceding section (n =1.7, 0) =0.2, 
p. =0.34738) the degree of destabilization is 20,/n =0.362." 

When 


2o,/n = 0, OF Mapp/n = 1 (31a) 
a system is completely stable. When 
0 <2¢,/n <1, or 1> tapp/n > 0 (31b) 
the system becomes partially destabilized (transition region), and when 
2c,/n>1, OF Mapp/n <0 (31c) 


the system is unstable. 

For conventional rotor blades m ranges from 1.4 to 1.9 (at sea level), and y~0.35 
represents the maximum attainable forward speed. By a systematic determination 
of o, associated with the values 1.4, 1.5,---, 2.4, 2.5 of m, and the values 
—1i, —0.9,---, +0.4, +0.5 of 6) the writer computed the degree of destabilization 
20,/n in ap, m region which covers the range of practical interest. The plots are shown 
in Figs. 3a, b. The calculated points are indicated by crosses. The points on the y-axis 
(small circles in Fig. 3b) were obtained by means of an auxiliary graph giving the plot 
of gand q’ vs. u for various n. In Fig. 4 this construction is shown for n =1.4, V3, 1.9, 
2.0, 2.1, 2.4. The intersections of the curves with hérizontal lines of ordinates 0 and 1 
yield the w values for which o, =0. 

A replot of Fig. 3, with uw as abscissa, m as ordinate and nqp,/n as the index of the 
curves leads to the stability diagram, Fig. 5. The diagram indicates two complete sta- 
bility regions, and two transition regions. The complete stability regions are charac- 
terized by app/m=1 and by response frequencies which vary with uw and m (between 
0 and 4w); the transition regions are characterized by the fixed frequencies 0 and }w, 
respectively, of the response, and by degrees of stability .,,/m which vary with p 
and m (between 1 and a minimum). 

From Fig. 5 it is clear that the flapping motion is very stable for 1 <0.5 at any n. 
This is in agreement with the conclusions reached in the preceding section for the spe- 
cial case of m=1.7. 





4 For n=1.7, 09=0 one finds » =0.65734, and 6; and 6. assume the values given in I, Sec. 5. The 
corresponding degree of destabilization is 20,/m =0.510. 










































































WILE 
ANNAN 
\\\ Se. 
\ Ae 
eal 

rf 


u/49ey -) ‘Noiwzimievis3a JO 334930 


(a) 























advance ratio, ua 



























































*NOLLYzI71¢ 





“ 
VLS30 JO 334930 


r=) 


_~ 
Ee) 
ed 


v2] 
& 
rr) 
> 
w 
= 
° 
‘= 
© 
> 
— 
i) 
° 
=) 
£ 
vv 
D 
a) 
= 
> 
ZU 
SS) 
w 
> 
a 
3 
wa) 
~ 
° 
2 
° 
a 
= 
° 
S 
° 
ee) 
NY) 
ons 
5 
5 
as 
ao) 
— 
° 
Y 
£ 
oa 
a) 
Q 
uv 
S 
Ss) 
As] 
* 
w 
2 
ba 


of the aerodynamic damping coefficient n. 
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Fic. 4. The functions g and q’. 
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Fic. 5. Rotor blade stability diagram. Level lines indicate the relative degree of stability for 
various values of the damping coefficient and the advance ratio. 
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It may be well to recapitulate at this point the various assumptions on the basis 
of which the stability of flapping motion was established. The assumptions are as 
follows: The blades are attached to the rotor shaft by means of flapping hinges only 
(the freedom of oscillation of the blades in the rotor plane, about the lag hinges, is 
thus disregarded), and the distance of the hinges from the rotor axis is negligible. 
The hinge axes are perpendicular to the rotor shaft and to the respective blade span 
axes. The blades are untwisted, and are rigid both in bending and in twist. The 
chord and the lift coefficient are constant along the blades. The uniform motion of 
the helicopter is not affected by the oscillations of the blades. Of the aerodynamic 
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Fic. 6. Stability chart of rotor blade flapping motion. 


forces which act on the blades only the velocity square lift forces are taken into ac- 
count. Air velocities in the direction of the blade span are neglected. The field of in- 
duced velocities is assumed as constant over the rotor disc, and interference effects 
between the various blade elements are disregarded. All angles involved (but W) are 
small, in particular the angle of attack is always less than the stalling angle. 

Although the above assumptions are violated by a helicopter in many respects, 
experience has shown that the steady state solution (10) gives a description of the 
flapping motion of the blades which, for 10.35, is in satisfactory agreement with the 
observations. Thus Eq. (5) can be considered as essentially correct. 

From the practical point of view the plots of Fig. 6, giving the percentage decrease 
per revolution of the amplitude of a transient, constitute the principal result of the 
present paper. On the one hand the chart establishes the hitherto unsuspected and 
curious behavior of the flapping transients near n =+/3 (the stability decreases as yu 
increases, the natural frequency of the flapping motion is 4w over a wide range of p’s 
and n’s); on the other hand it also focusses attention on the sluggishness of the blades 
in recovering their steady state motion when the value of n is small. Blades with low n 


must therefore be avoided.'® 


16 The widely held belief that “the natural frequency in flapping is equal to the speed of rotation” 
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The problem of blade response to control deflections leads to considerations which 
go beyond the question of the rate of extinction of a transient, and will therefore be 
taken up in a separate paper. 

5. The functional relation ¢,+i0;=o0(u, m). The characteristic exponent ¢ is a 
simple expression, (21), in the combination gq (or gq’) of 6) and D. @) and D in turn are 
functions of the advance ratio » and the damping coefficient m. As pw increases, 9, 
Eq. (19a), decreases monotonically. The variation of D with yw is more complicated. 
D becomes infinite for 69=0*, 1°, 2?, - - - ; for the flapping problem (@)<1) only for 
65=0. Elsewhere D fluctuates between positive maxima and negative minima.?* Cor- 
respondingly g maxima (¢, maxima) alternate with g’ maxima (¢, maxima). (Cf. 
Fig. 4 which shows the g(u), g’(u) variation near the origin.) Complete stability exists 
only in the very narrow strips where 0Sq3S1. In these regions the ratio map,/n is 
constant at the value 1, while the frequency varies across the region by }w. In the 
remainder of the uw, m diagram the dominant frequency of the transient is an integral 
multiple of 4w, and one has either 0<2¢,/n<1 (transition region), or 2¢,/n>1 (in- 
stability region). 

The relationships ¢,+10;=o(u, m) can be best appraised by considering the varia- 
tion of o with pw for one particular value of m. Such a survey is carried out below, first 
for n=2.4 in the w-range 0 to 0.5 with the aid of Fig. 5; then by an approximate 
method for n=+/3 in the u-range 0 to 8.0 wth the aid of Figs. 7, 8, 9. Although 
Eq. (5) ceases to give an approximate description of blade flapping when the stalled 
areas‘of the blades are large (when u=0.5), the equation is of interest, per se, also in 
this case, and an investigation of the behavior of the solutions of (6) for large values 
of u aids in the formation of a clear picture of the variation of o with uw and n. 

For n=2.4, u.=0, the transients 8; and 6: are solutions of Eq. (13), and thus can 
be written as 


8, = exp (— 34n + Vin? — 1)y = exp (— 1.2 + 0.663)y, (32a) 
Bs = exp (— 3n — Vin? — 1)p = exp (— 1.2 — 0.663)y. (32b) 


Therefore %app/n =0.447 and o;=0. The same result is obtained also from the general 
theory. Since for »=0 one has D=1, therefore 


q = sinh rV/ fn? — 1 = 3.95, (32c) 
1 oe 
o = — log (q + Vq’? + 1) = 0.663, (32d) 
T 
and again 


1 — 2¢,/n = 0.447, o; = 0. (32e) 


Tr 





(quotation from Bennett’s Princeton University Lecture Notes On the physical principles of the helicopter, 
p. 1, 1944) is largely the result of an unfortunate misuse of the expression “natural frequency” for the fre- 
quency at which the phase displacement between impressed force and response is 90°. For many slightly 
damped systems with constant parameters this frequency nearly coincides with the frequency of the 
transient,” but sometimes, as in the present case, there is considerable divergence between the two. 

17 J. P. Den Hartog, Mechanical vibrations, p. 63, McGraw-Hill, 1940. 

18 Oscillation theorem, cf. M. J. O. Strutt, Lamésche, Mathieusche und verwandte Funktionen in 
Physik und Technik, p. 15, Springer, 1933. 
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As p increases, g’ decreases until it reaches zero at u~0.334. Correspondingly 
a, also decreases to 0, the characteristic curve (the common boundary of the transition 
and the complete stability regions) is breached, g’ changes into g, and the complete 
stability region is entered. Between p~0.334 and w~0.352, g varies from 0 to 1. 
Correspondingly o¢,=0, while o; varies from 0 to 3. In the transition region that fol- 
lows o;=4, and ¢, varies from 0 to a maximum, then back to zero. (The diagram does 
not extend sufficiently to the right to show this variation.) 
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Fic. 7. Stability of the equation 8+nw8+u*(1+4nyu cos wt)8=0 by the “twin-ripple method” of 
approximation. Complete stability regions are blackened, transition regions (bounded by the curves 
Napp/n=1 and mtapp/n=0) are lightly shaded, and unstable regions are left blank. The encircled numbers 
indicate the dominant frequency of the transient (in units of w) in the respective instability and transition 
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regions. 


Fig. 7 shows qualitatively the conditions which hold for Eq. (6) at large values 
of uw. This diagram was obtained by replacing the variable coefficients of 8 and 6 in 
Eq. (6) by their constant averages first in the Y= —2/2 to +7/2 interval, then in 
the Y=7/2 to 3/2 interval. The resulting pair of differential equations was solved 
by the method of Meissner’® (better known as the Van der Pol-Strutt “rectangular 
ripple” method?*.?") which was adapted to damped systems by Shao Wen Yuan and 


the writer.?? 





19 Strutt, ibid., p. 39. 
20 Den Hartog, loc. cit., p. 389. 
21S, Timoshenko, Vibration problems in engineering, p. 164, D. Van Nostrand, 1937. 


2 G. Horvay and S. W. Yuan, Stability of rotor blade flapping motion when the hinges are tilted. Gen- 
eralization of the “rectangular ripple method” of solution, to be published soon. The procedures used in ob- 
taining diagrams like Figs. 7, 8, 9 are also discussed in this paper. 
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In Fig. 7 the regions of complete stability are blackened, the transition regions 
are lightly shaded, and the-regions of instability are left blank. Consider the ordinate 
n=*/3. For »=0 one has the solution 


Bi=cos}y, B2=sin}y, tapp/n = 1; (33) 


the transition region ® extends into the point (u, ”)=(0, ~/3) of the ordinate 
axis. As p increases, the stability decreases until a minimum of mep,/n~0.50 is reached 
around w~-0.85. Thereafter the stability again increases. 

In the present approximation of Eq. (6) complete stability (¢,=0) sets in at 
~1.24 and ceases at y~1.28. In the meantime o; has increased from } to 1. In the 
region o;= @) (the dominant frequency of tke transient is now w) o, varies from 0 
at u~1.28 to a maximum of ~1.0 at p~3.2 (yielding 1.,,/n~—0.16), and then back 
to zero at y~4.02. First comes the transition interval (1>.p»)/n>0) from w~1.28 
to u~2.27. This is followed by the instability interval (”.,,/n<0) from u~2.27 to 
p~3.83. A second transition interval (1>m.,,/n>0) from u~3.83 to up~4.02 com- 
pletes the (1) interval. There follows now a very narrow complete stability interval, 
p~4.02 to u~4.03, in which ¢, =0 while o; varies from 1 to 14. This is followed by 
the transition interval, this by the instability interval; and so on, indefi- 
nitely. Complete stability, transition, instability, transition, complete stability re- 
gions alternate, ad infinitum. 

The stability regions, associated with ¢, =0 and varying o;, become ever narrower, 
the instability regions, associated with o;=const. and varying o, become wider and 
graver. The orders ©) : ® @ , «+ Of the transition and instability regions are 
marked in the diagram. The regions 2), - + + can be called even ordered regions 
(the dominant frequency of the transient is 0, 2, 4,--- X}w); the regions 
@, 44), - ++ can be called odd-ordered regions (the dominant frequency of the 
transient is 1, 3, - - - X4w). The boundary curves m,»,/n=1 and 0 are also indicated 
in Fig. 7. Level curves which touch the 2 =+/3 horizontal line complete the picture. 

It is rather interesting to note that although Eq. (6) reduces to a damped Mathieu 
equation when uy is very small (u?~0), Fig. 7 shows a departure from the conventional 
stability diagram given for these equations. In the first place, the present diagram 
plots ” vs. u, whereas conventionally one plots @) vs. wu. More fundamental is the 
change in behavior near the origin. Since the oscillating terms of the spring force 
s(t)B tend to zero as n—0, no transition region, like the one at n=+/3, »=0, 00=}, 
wedges into the origin =u =0, 09=1. A second comment is also in order. The twin- 
ripple approximation used in obtaining Fig. 7 accounts in an approximate manner 
for the stiffness fluctuations produced by the $w*un cos wt-8 spring term of Eq. (6). 
In the process of averaging the coefficients of Eq. (6) both the wun sin wt-8 damping 
term and the @u2n sin 2wt-8 spring term are eliminated. However for large pu, the 
nu? term of s(t) has the principal effect. This effect can also be evaluated by a “twin- 
ripple” approximation, provided one first omits the two u-proportional terms of 7(8), 
and then averages the yp? term in the intervals 0 to 7/2 and 2/2 to x. The resulting 
stability conditions are indicated in Fig. 8 for the value n=+/3. : 

The diagram indicates that in the absence of the terms fwun sin yf and 
4w’un cos WB of F(8) the region of complete stability extends from 4 =0 to about 
u~1. For u>1 the system is mostly unstable. Comparison of Figs. 7 and 8 indicates 
that the high order instability regions crowd much closer to the origin than was sur- 
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Fic. 8. Instability troughs of the equation 8+nw8+w*(1-+-ny? sin 2wt)B=0, (n =+/ 3), by the twin- 
ripple method of approximation; napp/n =1 represents complete stability, 0<mapp <1 represents transi- 
tion, #app/m <0 represents instability. The calculated points are indicated by small circles. The encircled 
numbers indicate the “orders” of the regions. 
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Fic. 9. Instability troughs of Eq. (6) for n=v/3 by a “4-ripple approximation.” The calculated 
points are indicated by small circles. The encircled numbers indicate the “orders” of the regions. 


mised on the basis of Fig. 7. Moreover the instability caused by the yu? spring term is 
much graver than that caused by the first power p term. 
It is to be expected that in an exact solution of Eq. (6) the first power wu terms 
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determine the stability conditions for 4<1, whereas the yw? term is principally respon- 
sible for the stability conditions at y>1. This is born out by Fig. 9. Fig. 9 was obtained 
by use of a “4-ripple” approximation. (The interval 0 to 27 is divided into 4 equal 
parts and the coefficients p(t) of 8 and s(t) of B are replaced by their averages in each 
interval.) This approximation takes into account the effect of all the oscillating co- 
efficients of Eq. (6). The figure shows that for u>1 the odd ordered instability regions 
(which are produced solely by the y-proportional terms) are of minor importance as 
compared with the even ordered instability regions. 

6. Resonance.”* The preceding sections were concerned mainly with the transient 
solutions of Eq. (5). The present section deals with the steady state response. 

In general, the assumption (10) leads to an equation system (11) which can be 
solved for the Fourier coefficients ao, a1, - - - , Gn, by, - - -. When however the de- 
terminant A of the system vanishes, then the Ansatz (10) fails. In this case 


Bo = do + a; cos wt + 3; sin wi + a2 cos 2wt + - - - 
+ t[co + c: cos wt + d; sin wt + cz cos 2wt + --- | (34) 


is the correct form for the solution of the equation. Note that the response becomes 
infinite even though the system is damped. 

What are the conditions for A=0? Evidently, A=0 implies that the homogeneous 
equation (6) has a solution which is of the form (10). This in turn implies 


o=intmi (m=0,i1,2,---). (35) 


Thus, resonance occurs on the napp/n=0 boundary lines of the even ordered transition 
and instability regions. From Fig. 5 it is seen that for helicopters resonance excitation 
is just as remote as instability.™ 

There are two more interesting questions which arise in connection with resonance. 
(1) How would conditions change if the exciting force on the right side of Eq. (5) were 


of the form 


sin 
E(t) = { bn, where v * w, 2w,--- ? (36) 
cos 


(2) How would conditions change if the system were undamped, i.e., if the p(é)8 term 
were absent from Eq. (5), but the s(#)8 expression still contained oscillating terms of 
frequency w? 

If y=4w or an odd multiple of it, then placement of 


Bo = a; cos 4wt + 5; sin Fwt + a3 cos Fut + 53 sin foi + --- (10’) 


23 This section is based on the paper by G. Kotowski, Liésungen der inhomogenen Mathieuschen Differ- 
entialgleichung mit periodischer Stérfunktion beliebiger Frequenz, Z. angew. Math. Mech. 25, 213, 1943. 
Kotowski credits introduction of the method to A. Erdelyi, Arch. Elektrotechn. 29, 473, 1935. 

* Since the response amplitudes are large near a resonance frequency the following question is of 
interest in many problems: how deep do the regions of large amplifications extend into the transition re- 
gions or perhaps stability regions? Kotowski answers the question for the equation y’+(A-+cos x)y 
=a4-a; cosx+az cos 2x, by plotting the resonance curve. (A is the variable parameter.) In the present 
problem the great distance of the helicopter operating range, u<0.35, from the nearest napp/m =0 curve 
eliminates the necessity for such an investigation. 
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into ¥(8) = E(t) leads to an equation system (11’), with a determinant A’ which in 
general does not vanish; A’ vanishes only when the homogeneous Eq. (6) has a solu- 


tion of the form (10’), i.e., when 
o=in+ (4+ m)i (m = 0, + 1, +2,--:) (35’) 
Thus, for v=odd multiple of 4w, resonance, occurs on the napp/n=0 boundary lines of 
the odd ordered transition and instability regions. In this case expression (10’) plus ¢ 
times a similar expression is the correct ‘“‘Ansatz”’ for the forced response. 
When » is a fractional multiple of }w, then 
Bo = cos vt[ap + a; cos wt + b; sin wt + az cos 2wt + --- | 
+ sin vt[co + cy cos wt + d; sin wt + ce cos 2wt+--- | (10’’) 
is the correct assumption for the particular integral. The associated equation system 
(11’’) has a determinant A’’ which vanishes only when the homogeneous equation (6) 
has a solution of the form (10’’), i.e., when 


Vv 
c= int (—+m)i (m=0,+1,+2,---). (35”) 
w 


But from Eq. (21) it is obvious that o is never of this form, unless n =0. Thus, when 
2v is a fractional multiple of the parametric frequency w, resonance can occur only when 
the system is undamped. In such a case the response occurs along the lines o,=m+v/w 
of the complete stability regions, and (10’’) plus ¢ times an expression of type (10’’) is 
the correct form for the solution. 

It is to be noted that a single frequency v excites an infinity of response frequencies 
v, vitw, v+2w, - +--+, and thus in a system with variable spring and damping coeff- 
cients one obtains an infinity of resonance peaks in contrast with a system character- 
ized by a 2nd order differential equation with constant coefficients which shows only 


one resonance peak. 
The absence of damping also modifies the results for the cases when rv is an even 


or odd multiple of 4w. The 2.,,/m=0 curves now coincide with the characteristic 
curves of the stability diagram, and 


D> [Cae + bat + cxf*) cos kat + (de + ext + fet?) sin kot] (37a) 
0 

and 
> [Caz + bet + cut?) cos (2k + 1)ot + (de + ext + fat?) sin (2k + 1)wt, (37b) 
0 


respectively, are the proper forms for the resonant response. 
The proofs of the expansions (34), etc. are based on the well-known formula” 


Y . 
Je"(W)Bo(y) = uv) f v3(x)e (x) E(x)dx — uv) f vo(x)e“"(x)E(x)dx (38a) 





% There is no need to write out the equation system. 
% Frank-Mises, Differential und Integralgleichungen der Physik, Vol. 1, p. 300, Fr. Vieweg u. Sohn, 


1935. 
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where e(W) is given by Eq. (17a); 2:(W), v2(W) are the two linearly independent solutions 
of Eq. (18); 
J = 2,07 — vev{ = const.; (38b) 


and E(y) is the exciting function. 
When 2, is of the form 


e-"¥/2 > geikv (39a) 
and the exciting function is of the form 
- cos) . : 


then the integrand 
v2(x)e“"(x) E(x) 
is a series which contains a constant term. This term integrates into y, and, multiplied 
by —e(W)v:(W)/J, yields the non-periodic 
t>> [cx cos ky + dy sin ky] 
terms of the expansion (34). The remainder of (38a) yields the periodic terms of (34). 


This proves (34). The proof of the other expansions is similar.?’ 
When the exciting function is of the form 


Ey) = ¥ cos (40) 


then, by (38), the response is of the form 
Bo(W) = D> [(ax + cap) cos ky + (bx + dip) sin ky] (41) 
0 


Substitution of (41) into the non-homogeneous equation (8) = E(t), and a compari- 
son of the coefficients of cos ky, sin ky, ¥ cos kw, y sin ky on right and left, yields 
an infinite system of equations (11’’’) for the coefficients a,, bx, cx, dx which can be 
solved in the manner of Eqs. (11). Use of the expansion (41) will be made in the study 
mentioned at the end of Sec. 4. 





27 The details (for Mathieu’s equation) are given in Kotowski’s paper. Kotowski also points out that 
there are certain exceptional cases. For instance, when a homogeneous Hill equation has a solution of the 
form > C; cos ky but none of the form > C; sin ky, then the latter expression is admissible as forced re- 
sponse to an excitation E(y) =sin y. 








168 


TORSION OF A CIRCULAR TUBE WITH 
LONGITUDINAL CIRCULAR HOLES* 


BY 
CHIH-BING LING 
California Institute of Technology 


Introduction. This paper presents a solution of St. Venant’s torsion problem of a 
circular tube having a ring of uniformly distributed longitudinal circular holes of equal 
radii. The investigation of the stress distribution in such a tube is not only of theo- 
retical interest but also of practical importance. The solution is found by constructing 
three series-of harmonic functions which are invariant with respect to rotation about 
the axis of tube. The parametric coefficients attached to these functions are deter- 
mined so as to satisfy the required boundary conditions. The case in which the central 
hole does not exist is also investigated. Finally, the solution is illustrated by working 
out several numerical examples. 

It may be mentioned that the tor- 
sion problem corresponding to the case 
without a central hole was solved by 
M. Kondo a few years ago by using the 
so-called stress distribution method.' 
The same method has been used by 
Th. v. Kérmd4n with success to find the 
potential flow around an airship, but 
its application to the torsion problem 
by Kondo seems to result in great com- 
plication. In Kondo’s paper, no formula 
for the torsional stiffness is given. 

Method of solution. Consider a cir- 
cular tube defined by the exterior and 
the interior radii @ and c respectively. 
Let there be a ring of longitudinal cir- 
cular holes parallel to axis of the tube, 

Fic. 1. Cross section of tube. each of radius bd with their centers uni- 
formly distributed over the circumfer- 
ence of a concentric circle of radius } at points represented in the Argand plane by 


z= be? milk (1) 


Y 





where 3 is the distance from the center of the tube to the center of any eccentric hole, 
k is the number of eccentric holes and m=0, 1, 2, - - - (@—1); the origin being at the, 
center of the tube. Fig. 1 shows the cross section of the tube. 

Define the following two systems of polar coordinates 


z= x+ iy = are‘, be = z — b = Dpe*, (2) 
where ar, ¢ are polar coordinates referred to the center of the tube as origin and bp, 6 


* Received Sept. 18, 1946. 
1 Masaiti Kondo, The stresses in a twisted circular cylinder having circular holes. Phil. Mag. (7) 22, 


1089-1108 (1936). 
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are polar coordinates referred to the center of one of the eccentric holes as origin, 
r, p being dimensionless. Also, for brevity, denote the ratios of the radii by 


b Cc 
pees: See (q<p< 1) (3) 


a 


where a is the radius of the exterior boundary of the tube and c is the radius of the 


central hole. 
It is observed that the solution of St. Venant’s torsion problem of the tube in 


question requires a harmonic function y, with its conjugate function single-valued, 
which satisfies the following three boundary conditions: 
(1) on the exterior boundary of the tube where |z| =a 


Y= y— (x*+ y9)/2 = 0; (4) 
(2) on the boundary of the central hole where | z| =¢ 
WY = const. = a; (5) 
(3) on the boundaries of the ring of eccentric holes p=A 
Vv = const. = a, (6) 


where bd is the radius of the eccentric hole under consideration, \ being a dimension- 
less quantity. Besides, by symmetry, the function y must be invariant with respect 
to rotation about the origin through angles 2m7/k and must be even with respect to 
the lines of symmetry. 

Hence we write 


Y=entvbet+Vs (7) 
and construct 
vi = a?) A,r" cos nko, 
n=0 


oe: 


v2 = a?), ; COs nko, (8) 
n=1 ie 

vs =a?) C.U,, 
a=] 


where the quantities A, B and C are dimensionless parametric coefficients to be ad- 
justed so as to satisfy the required boundary conditions. The class of harmonic func- 
tions U,, first obtained by R. C. J. Howland? is defined by 
Wo = Uo = iVo _— = log (s* — b*), 
b* d*Wo (9) 
(s—1)! db* 





W, = U, — iV, 


A brief account of these functions will be found in Appendix A. Note that in the last 


2 R. C. J. Howland, Potential functions with periodicity in one coordinate, Proc. Camb. Phil. Soc. 30, 
315-326 (1934). 
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Eq. (8) the function U> is rejected on the ground that it gives rise to a multi-valued 
conjugate function. 

Now, to apply the first and the second boundary conditions, we first express ws 
in terms of the polar coordinates (r, @) referring to the center of the tube as origin. 
The expansion of W» into a power series is different according as | z| is greater or 


smaller than 6. When | z| >b, 
- 1/6\* 
- blogs +X —(—) , 


Wo = 
n=1 Zz 
and when | z| <b, 
ead 1 Z nk 
Wo = ri— klogb+ >> —(=) ‘ (10) 
ant #\O 


Differentiation gives 


~ f(nk — 1 b\"* 
w.=k D(’ )(=) , 

aks] z 

“ (nk+s—-1 s\n 
— $)°s — 
(“4 x ( s—1 \) . 


respectively. It is noted that the binomial coefficient (}) vanishes when s>n. Taking 
the real part and observing that p=b/a, we have 


~ (nk —1 nk 
U, k>, ( \(4) cos nkd, 
n=1 | 1 rT 


(11) 
W, 





~ f/nk+s—-—1 r \nk as 
U,=(- ne d( )(<) cos nkd. 
n=0 - = 1 p 
Hence the last Eq. (8) becomes 
rs Pe | 1 ) 
vs = a’), —~ cos nkp, when |2z| >) 
cae (13) 
vs = a? > Qnr™* cos nko, when | z|<b 
n=0 
where 
nk (nk — 1 
p= tt B(P~ 
e=1 » Sa 1 
(14) 
k > 1) - 
Cn i pn* s=1 s-—1 “i 


Applying the boundary condition (4) at |z| =a or r=1 and using the first series 


for W3, we find 


1 eo 
O= 44—-— + > (An + Ba + Px) cos nko. 


n=1 
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This relation is satisfied, provided that the coefficient of each term vanishes identi- 
cally, i.e., that 


=0, A,t+tBitP,=0. (15) 


vl 


Ap — 


Applying the boundary condition (5) at | z| =c or r=q and using the second series 
for ¥3, we obtain 


Wi = Ao + Qo — g?/2 + DE (Ang™* + Bn/g™* + Qng”*) cos nko. 
n=l 


Equating the constants from both sides and setting each coefficient equal to zero, we 


find 
Ap + Qo — g?/2 = VW, Ang™* + B,/q"* + Qrg™* = 0. (16) 


To apply the third boundary condition, we first transform ¥V in terms of the polar 
coordinates (p, @) referring to the center at (x, y) =(0, 0) as origin. Thus, 


x? + y? = 57(1 + p* + 2p cos 8) 


and 


nk /nk 
r™* cos nk = R(z/a)"* = p**R(1 + pe®)"* = p*>> om cos mé, 
m 





m=0 
cos nko a\** 1 me of 
—— = Ri — =. = R(1 + pei?)—nk (17) 
rk Z ae 
. = nk-+m-— 1 
= } (— 1)” p™ cos mé. 
i m 


Consequently, the first two Eq. (8) become 


vi = a? >> Mnp™ cos m0, = 2 = 2.) Nunp™ cos md, (18) 
m=0 


m=0 


where 


3 k = k —1 
es o("") prt Na@ (- y= (" vi ) 8/9" (19) 


n=0 7 r=1 


Again, according to Appendix A, Eq. (A.12), 


U, = cos s0/p* + ps (— 1)" "a,p" cos n6. 


Thus, 
how {be + D (Calo + Imp”) cos mo | (20) 

where 
ite thin > mo,Co. (21) 


s=1 








CHIH-BING LING 


Applying the boundary condition (6) at p=A, we obtain 


=M,+N+h- p(1 a d?) /2 _— pr cos 8 


+ > (MadA™ + Nad™ + Cn/X™ + Lnd™) cos md. 


m=1 
Equating the constants and setting each coefficient equal to zero, we find 
Mo + Not Lo — pl + dr?) /2 Vo, 
MA+ NA+ C/A + LA — pA = 0, (22) 
Mar* + Nad* + Ca/dA” + Led”* = 0. 


From the first Eqs. (15), (16) and (22), we see that 
Ao = 1/2, VY; = Qo + (1 = q*)/2, We = Mo + No + La —_ Pr + r*)/2 (23) 


WV, and ¥, are thus determined, if the parametric coefficients can be found. The sys- 
tem of simultaneous linear equations represented by the other equations is just sufh- 
cient for the determination of all parametric coefficients involved. They may be re- 
written as follows: 


as + B, + | age = 0, Ae + pions, + — 0, 
q 


(24) 
A~*C 1 + M, + N, - Ly = p?, 173°C + es + Na + , i = 0 
The first two equations give, for n21, 
i= (q?"*O, _ P,)/(1 i g?"*), B= g?"*(P, cee Q,)/(1 q?"*), (25) 


Substituting these values into M, and N, in the last two equations, we find with the 
aid of (19), (14) and (21) 


= pr? + >> BL, Cu=d}9"> BC. (26) 


s=1 e=1 
where, in general, 


™B, = — (-— 1) m mo. 
sel yi" - ' ( q Ae +s- | 
k {- ane at gs 
7 > g?"* \m s—1 ‘ p es 
(— 1)"q?** /nk+m-—1 e i= ') , Ee = ‘yf 
ae a. eos 1 P — 2nk . 27 
: i oe ne m ik s—1 4 ent EB si 


This system of linear equations can be solved by successive approximations as follows. 
Write 


jcc. (28) 


where 
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and, by iteration, 


Ce =r" > B,C, ”. (29) 
e=1 


With the values of C,, thus found, the coefficients A, and B, can be obtained without 
difficulty from (25) with the aid of (14). 

It is to be observed that the validity of the solution by successive approximations 
depends upon the convergence of the series (28). From physical considerations alone, 
it seems likely that there will be convergence as long as the boundaries do not over- 
lap, i.e., 

q+ pr<p<1-— pr, A < sin (4/k). (30) 
To establish convergence of such a series analytically is usually difficult. However, 
we may proceed by first using some inequalities for the coefficients so that the series 
(29) can be summed, and then applying the ratio test to ascertain the range of con- 
vergence. Obviously some convergence is lost in thus using the inequalities, but it is 
often possible to establish a considerable range of convergence in this way. In the 
interest of brevity, no details will be given here. 


Torsional stiffness. The torsional stiffness of the tube is given by 
H= 2f f Wdxdy + 2nc?a’¥, + 2krd*b?aWe, (31) 


where the double integration is extended over the entire cross-sectional area of the 
tube, excluding the holes. Let 


H;= 2f yjdxdy, (j = 1, 2, 3), HW, -f (x? + y*)dxdy. 
Then, 
2f [ vaxdy = H,+H8.:+4H;3-— H,. (32) 


These integrals, save H3, can be evaluated by ordinary methods of integration with- 
out difficulty. The following results are obtained: 


H, = ra‘(1 — g? — kp’r*?), He = 0, (33) 


Hy = wat{1 — qt — kp*n(2 + d2)}/2. 


The integration of H; by ordinary methods fails since there exist singularities at 
the centers of the eccentric holes. However, H3 can be evaluated by using the follow- 
ing relation derived from Green’s theorem, in which F is a function of the complex 
variable z regular in the domain S enclosed by a contour C. 


dF 1 
ff — dS = —if Fdz, (34) 
Ss dz 2 “4 


the contour C being taken in a counter-clockwise direction. In this equation set 
W,=dF/dz, so that 








174 CHIH-BING LING [Vol. V, No. 2 


F -f W,dz, (35) 


and let the contour C be taken around the sector of 1/2k of the cross section, as shown 


in Fig. 2. Then 
ff U,dxdy = 2er f f W,dS = Rif Fd3. (36) 
Ss c 


Using the appropriate series of W, in 
different parts of the contour, we find 
that this integral vanishes identically 
for all values of s except for s=1. The 
term arising from 1/f¢ in W, does not 
vanish when the integration is taken a 


along the semicircle p=X. Here - 
4 \ 
at 4\_\ 
0 


F fs biog f 
= —_—_ = og P 
¢ 


and when p=A, 






Fic..2. Contour C of 1/2k of section. 


F = d(log \ + 7), dz = bdde~*®, 


We thus find 


0 
f Fdz = b*\(2 + 2 log A + ix). 


Hence, 


ff U,dxdy = — rkyb?, (37) 


H; = 2a? >) cf f U,dxdy = — 2wkdb?aCy. (38) 
s=1 


and consequently, 


Summing up, we find the torsional stiffness 
H = ra‘{(1 — g°)? + ptr? (2 + d*) — 2kpPr (2C1 + A) + 4g + 4kp’NWe}/2. (39) 
The resulting twisting couple is given by 

T = pH, (40) 


where yu is the modulus of rigidity of the material and 7 is the angle of twist per unit 


length of the tube. 
Stress components. The non-vanishing stress components of the tube are two 


shear stress components given in terms of the polar coordinates (r, ¢) by the following 


° = — «——, g eS a= es opem's (41) 
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while in terms of the polar coordinates (p, @), they are 


ov ov 
eee eee i ide Neat ee (42) 


— , 26 
bp 06 b dp 


It is now a straightforward matter to find the stress at any point of the cross section. 
By using appropriate series, the shear stresses on the three boundaries are found as 





follows: 
[25 ]n1 = pra (1 — 2k > nA, COS nko), 
n=l 
2 > we 

[26 ]p=n = pra (mr + Dr » 7 cos nt), (43) 
2 <=. af. 

[26 ],no = mura (« +— >> —— cos nko). 
q naa’ qv 

fs OO 
CASE I CASE IL CASE I 


Fic. 3. Illustrative examples. 


The maximum shear stress occurs on the boundary across the narrowest sections. 
Due to the presence of the holes, the greatest shear stress does not necessarily occur 
on the exterior boundary. 

The case without central hole. It is interesting to investigate also the case in which 
the central hole of the tube does not exist. A solution is obtained for y by merely 


omitting Y2 in Eq. (7), such that 


¥ = + Ys. (44) 
Consequently, the relations obtained from the boundary conditions become 
Ao = 3, WV. = Mo a Lo —_ $ 2(1 oa A*), (45) 
A, = — P,, 77C,+ M1+ l= fp’, 12C".. + Mp +L. = 0, 


so that the coefficient in (27) is simplified to 


= k k—1 
mg, = —(—1)™ "a, + kX (" \(" yom, (46) 


n=l \M s—1 
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Numerical examples. The foregoing method of solution will be illustrated by work- 
ing out the following three examples: 











Case k p q a No. of Holes 
I 2 } 0 3 2 
I] 2 3 } 3 3 
III 6 3 } 3 7 


The arrangements are shown in Fig. 3. In each case the holes have equal radii and the 
nearest distance between any two adjacent boundaries is equal to the radius of one 
hole. These cases probably represent the most crucial ones ever met in practice. 
Tables of the coefficients "8, are first prepared, with which the parametric coefficients 
C, are found by the method of successive approximations; the results being shown in 
Appendix B and in the accompanying Table I. The shear stresses along the bound- 
aries are then computed and shown in Table II. Fig. 4 shows the results diagram- 





























matically. 
; TABLE I. The coefficients C, 

m Case I Case II Case III 

1 7.8560 X 1072 4.2392 K 107? 3.0593 KX 107 

2 3.1549 x 10-3 3.1240 10-4 5.8718 X10 

3 —3.0714x10- } 1.9115 10-5 4.9273 1075 

4 9.6791 K10-5 3.2509 X 1077 —7.4587 X 10-8 

5 —1.7036X107 1.8686 1077 4.9622 x10-7¢ 

6 | 3.9811 10-8 —1.2587 10-8 5.9963 x 10-8 

7 —8.0836 1077 2.5340 x 10-9 — 1.0635 «10-8 

8 1.7422 x1077 —2.7810 «10- 4.7674 10-1 

TABLE II. Shear stress/ura along boundaries. 
| Case I Case I] | Case III 
ke |__ eae te ee ee ee) ee ce 
Or | Exterior | Eccentric | Exterior | Eccentric | Central Exterior Eccentric | Central 
° boundary holes boundary holes hole | boundary holes hole 
0 1.185 Lp |. 3 308 1.547 0.459 15327 1.247 0.214 

30° 1.13 1.049 1.124 1.327 0.402 | 1.105 1.083 | 0.212 
60° 1.022 0.660 | 0.967 0.793 | 0.277 1.050 0.583 0.207 
90° 0.956 0.198 | 0.930 0.145 0.160 | 0.988 0.068 | 0.200 
120° 0.928 —0.199 0.928 —0.451 | 0.081 0.938 —0.259 0.193 
150° 0.918 —0.396 0.931 —0.871 0.039 0.907 —0.529 0.188 
180° 0.916 —0.392 0.932 —1.033 0.026 0.897 —0.678 0.186 


It appears that in all the cases the greatest shear stress occurs at those points on 
the boundaries of the holes which are nearest to the exterior boundary. The magni- 
tudes are equal to 1.194, 1.547 and 1.247 respectively except for a factor ura. The 
points of greatest shear stress will, however, shift to the exterior boundary if the 
eccentric holes are sufficiently rer ste from the exterior boundary. It is seen that in 
each case there exist two null points or points of zero stress on the boundary of each 
hole at an angle @ numerically between 90° and 120° or, to be more precise, at 105°, 
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96° and 95° respectively. It is also seen that due to the presence of holes the maximum 
shear stress on the exterior boundary is increased by 18.5%, 30.8% and 12.7% re- 
spectively. From the stress distribution on the boundaries it appears that among 
the three cases the case II is the most unfavorable arrangement in which, it is re- 
called, the centers of the holes are situated collinearly on the same diameter of the 
section. The stress distribution on the boundary of the central hole is less important 
since the magnitude is smaller. 

Table III shows the cross-sectional area A and the torsional stiffness H of the tube, 
in which the last two columns indicate their percentage reductions due to the presence 
of holes respectively. The relative reductions in torsion stiffness due to 1% reduction 
in area are equal to 0.729%, 1.589% and 0.963% respectively. Again, it appears that 
the case II is the most unfavorable arrangement among the three where the effective- 
ness of torsion modulus is concerned. 


TABLE III. Cross-sectional area A and torsional stiffness H. 














| 
Case A/rxa?* H/}xa* % reduction of A | % reduction of H 
ae © } | 
I 0.8367 | 0.8809 16.33 11.91 
II 0.8800 | 0.8093 12.00 19.07 
Il 28.00 26.95 
{ 








The author wishes to express his deep appreciation to Dr. H. S. Tsien of the 
Guggenheim Aeronautical Laboratory, California Institute of Technology, for helpful 
criticism and suggestions. 

APPENDIX A 

Class of functions invariant with respect to rotation. Let there be on the circle 

|z| =b in the Argand plane a series of k points 
ga §ert?, (A.1) 


where m=0, 1,2, +++, (R—1). 
A function with a logarithmic singularity at each such point is defined by 


k-1 
Wo = — log (2* — b*) = — log (2 — b) — D> log (2 — be2™*4/*), (A.2) 
m=1 
Writing z—b=d¢, we have except for a constant 
k—-1 
Wo = — logt¢ — >> log (1 — umf), (A. 3) 


m=1 


where 
1/tn = 1 — e®eril’, 


such that u» is a root of the equation 
(u— 1)* —u* =0. (A.4) 


Expansion yields 


Wo = — loge + > (— 1)* "ang", 
n= 1 








eh ee 
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where 


1 k—-1 m 
"a9 = — > un = o,/N, 
nN m=1 
¢, denoting the sum of the nth power of the roots of Eq. (A.4), which may be calcu- 
lated for any given values of k and n. 

It is readily shown that functions derived from Eq. (A.2) by differentiation with 
respect to b-are invariant with respect to rotation about the origin through angles 
2mr/k. If we define 

: b* d*Wo e 
W,= , (A.7) 
(s—1)! db* 





we find 
1 1 (1— 1)" 1 k-1 - n+s—1 
W,=—+ —— =— + tn — 1)° 2, (- p*( ) tut)”, 
, ~ a+)? fF pA n ( 


$ 





or 


1 oe 
W, = a + >) (= 1)" "ag", (A.8) 


n=0 


where 


"a, = (" sy p> (tim — 1)*tm = (" = ') Fi. D*(* Jones 
n 


n m=1 t=0 
n+s-—1 
= ( ) atom (A.9) 
n 
A’o, being the sth finite difference of the series On, On41, On42, °* *- From Eq. (A.4) 


we obtain 
mk mk 
yrtmke — = (-— 1)! ( ) yrtmk—t 
t=0 t 
where m is an integer. Summing over the roots, 
mk mk 
onimk = > (- y*( ) cnt = A™*o,,. 
t=0 t 
Hence 
Am*teg, = A*A™*c, = A*Gaimt; (A.10) 


so that all the coefficients can be found from the values of ¢, and its first (k—1)th 
finite differences. 
Now, if we define 


Wo = Uo — iVo, W, = U,— iV, (A.11) 


and ¢ =pe*, we find by equating the real and imaginary parts 





Il 





| 


ne: ie ace 
a 


ob ob ob 





—1.4995 x10 


—1.1196x107 


—1.0842 «107? 


ONDA MP Wh = | 


7) 


—9.9512 X107} 


| OIAMNR WHE | 


| 
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4 


= — logp+ p (— 1)” "aop" cos nO 


n=1 


g— p> (— 1)" “aop” sin 8, 


n=1 


i) 
zz (— 1)” "a,p" cos né 
n=0 
co) 

= >> (— 1)" "a,p" sin né. 
n=1 


APPENDIX B. TABLES OF COEFFICIENTS al 
Case I. k=2, p=3/7, g=0 








2.0736 107 —5.1703 X10 
7.2751 X 107 1.6111 «107 
2.4167 X10 —1.1142«107 
1.2817.x1071 1.8679 X107 
1.2605 X107 —1.2240 «107 
7.8758 X 107? 1.1563 107! 
5.5633 K 107? —7.7962 X107? 
3.4860 x 107 5.9407 x 107? 
m=6 m=7 
7.9832 x 10-3 —3.8723 X1073 
2.6253 X 107? 1.5891 K10~? 
5.7816 X 107 —3.4184 «107 
—7.8946 «107? 5.5966 X 107 
1.0893 X 10-1 —7.7443 X10 
—1.0718 10" 1.0010 «107 
1.1678 X10— —1.018010— 
—1.0237 X10" 1.0703 X10- 
Case II. k=2, p=3/5, gq=1/5 
m=2 m=: 
5.8617 X107 3.2014 K107 
1.3668 9.7830 X 107! 
1.4633 2.9506 
2.4953 2.3101 
4.4025 107 5.6549 
3.1571 9.7328 xX 107 
—1.6578 7.6527 
4.4190 —2.9699 





= ome w 


_ 


Noe F&O RK BS 


_ 
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(A. 12) 


m=4 


.4879 X10 

.4085 x 107? 
.4009 x 10-4 
.1476 107! 
.517010 
.1495 x10 
.0623 x 107 
.9025 K 10-2 


m=8 


.9615 X 10-3 
.7149 X10-% 
.2278 X 107 


.9978 X 107 


.1346 X 1077 


. 7385 X 107 


.3473 X10 
.6648 X 107? 


m=4 


.0272 X10™ 
.2483 
.7326 
.1916 
.1153 

.2122 X10 
.5465 
.7278 X10 
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$s m=5 | m=6 m=7 | m=8 
1 2.5737K10, | —1.6585x10-2 | 2.8770X107 | —2.9011x107 
2 1.761110 | 1.0524 | —4.6363x107 | 1.1047 
3 3.3929 | 4.866410? | 3.3914 | 4.1137 
4 3.2923 | 8.0811 1.4551 8.6392 
5 1.291710 6.6271 1.8327x10 | 4.0214 
6 7.9600 | 2.6963 X10 9.7852 4.073410 
7 2.521210 —s| 1.6097X10 | 5.659510 | 2.917810 
8 7.3042 | 7.431210 | 2.6870X10 | 1.196210 
Case III. k=6, p=3/5, g=1/5 
s m=1 m=2 | m=3 | m=4 
2.8379 | 1.6559 | 1.0945 | —1.7600 
2 3.3118 | —5.8522 | 5.9408 | 5.8842 
3 3.2836 | 8.9112 | —1.7916x10 | 1.839210 
n —7.0400 1.1769x10 | 2.452310 —6.6469 10 
5 5.6284 —2.8311X10 | 3.8057X10—s| 7.6789 X10 
6 6.0363 2.2409 x10 | —1.087110® | 1.3393 X10? 
7 —1.3620x10 | 2.8477 X10 8.8391X10 | —4.1222x102 
8 7.9354 | —7.0520x10 1.214810? | 3.4151 X10? 
5 m=5 | m =6 m=7 m=8 
1.1569 1.0060 —1.9457 | 9.9193 X107 
2 —1.132410 7.4698 8.1366 | —1.7630X10 
3 2.283410 —5.4335K10 | 3.7881X10 | 4.5598 X10 
4 6.143210 8.9290 x 10 | —2.3555X10? | 1.7075 X10* 
5 —2.4343 X10? 2.215010? | 3.4389X102 | —9.8099%102 
6 2.658110? | —9.0308x10® | 8.1328X10? | 1.3199 108 
7 4.8145 X10? 9.4883x102 | —3.3856x103> | 3.0419 108 
1.7598 x 108 3.4764 108 | —1.2761 104 


8 —1.5696 X 108 
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ANALYSIS OF THE TURBULENT BOUNDARY LAYER FOR 
ADVERSE PRESSURE GRADIENTS INVOLVING 
SEPARATION* 


BY 


W. S. COLEMAN 
Research Department, Blackburn Aircraft 


Part 1. 


1. Introduction. In a recent paper by von Doenhoff and Tetervin,! and subse- 
quently in another by Garner,? a new differential equation of an empirical character 
has been advanced, whereby the analysis of turbulent boundary layer development, 
and the prediction of turbulent separation, are greatly facilitated. The present paper 
deals with the problem theoretically, and, from considerations of the fundamental 
equations of motion, establishes an analogous expression which, though analytically 
different from that proposed by the above authors, contains the same parameters. 
It is further shown how this equation may be solved numerically by an approximate 
method of a rather complex character. 

On comparing results with those obtained from the empirical relations, it appears 
that there is good agreement when the pressure gradient is small, but that discrepan- 
cies are more serious for larger gradients tending to separation. In this respect there is 
a notable divergence between the results given by von Doenhoff’s formula on the one 
hand, and that of Garner on the other, the former being appreciably smaller at values 
of H in the region of separation. It appears, however, that Garner’s relation involves 
the three basic parameters of the theoretical equation, whereas, in the case of von 
Doenhoff and Tetervin, only two are apparent. There is reason to think, therefore, 
that Garner’s treatment may be the more reliable, a conclusion which is supported 
by the theoretical calculations, in that they are in very much better agreement with 
Garner’s predictions than with those of von Doenhoff and Tetervin. Nevertheless, 
despite the above discrepancies, step-by-step integration of the empirical equations, 
in conjunction with the momentum equation, leads to estimates of the boundary layer 
characteristics in -good agreement with experiment, so that the errors appear to be 
less important in the final result. 

So far as the theoretical treatment is concerned, it would seem of value in estab- 
lishing the essential parameters associated with the new equation, but for ease and 
rapidity of calculation in practical cases the empirical approach may prove more 
attractive. 

We shall now proceed to elaborate the fundamental arguments leading to the 
equation concerned. 


* Received March 22, 1946. 

1 A, E. von Doenhoff and N. Tetervin, Determination of general relations for the behaviour of turbulent 
boundary layers, Nat. Ad. Comm. for Aeron. confidential report No. 3G13 (1943). Also reprinted by Aeron. 
Res. Comm., Report 6845, F.M.597, Ae. 2255 (1943). 

2H. C. Garner, The development of turbulent boundary layers, Aeron. Res. Comm., Report 7814, 
F.M. 705 (1944). 
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2. Equations of fully developed, turbulent flow. The Reynolds equations? of tur- 
bulent motion in a two dimensional field of incompressible, viscous flow take the 
well known form 


Di 1fa stadies 
Di -—|=G.. - puu’) wae = (Bux — 9 u |, (2.1) 
Di 1fa =—s atin 
Di = a: E (Diy — pu'v’) + = (Puy — aa) | (2.2) 


There is also the equation of continuity, 
Ou 80d 
af nt (2.3) 
Ox oy 
where bars denote temporal mean values and primes the fluctuating components. 
Adopting Boussinesq’s procedure‘ of representing the Reynolds stresses pu’u’, 
pu'v’, pu'v’ by an apparent increase of viscosity, we replace the natural coefficient 
of viscosity u by (u+pe), where ¢€ is the measure of the momentum interchange due to 
turbulence, and must therefore be considered to vary with respect to the space co- 
ordinates. The Reynolds stresses are further treated in exactly the manner prescribed 
by the general theory of stress. Accordingly, 





—_ Ou 
pez bod pu'u’ va eee p si 2(u > pe) yO 
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whence the bars denoting mean quantities are no longer necessary, and will therefore 
be omitted in what follows. 

3. Curved flow. Intrinsic form of the equations. Equations (2.5), (2.6) will now be 
referred to the curvilinear co-ordinates of Fig. 1. Let AB be a segment of any stream- 
line of the mean flow, and let PN be the orthogonal curve through the point P on the 
streamline where the resultant velocity is g. Let ds, dn be elements of arc of AB and 
PN respectively, and Rneny. let @ be the angle the tangent to PN at P makes with the 
axis of X. 

20. Reynolds, On the dynamical theory of incompressible viscuous fluids and the determination of the 


criterion, Phil. Trans. (A) 186, 123-164 (1895). 
‘T. V. Boussinesq, Essai sur la théorie des eaux courantes, Mém. Sav. Etrang. 23, No. 1 (1877). 
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and the equations of motion for steady flow take the form 
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4. Application to the boundary layer. Integral equation for 6, In applying the equa- 
tions of motion to the boundary layer, we make the usual assumptions that 

(a) the fluid is of small viscosity; . 

(b) the boundary layer thickness 6 is small, O(v'/) ;* 

(c) the curvatures 00/0s, 00/dn are nowhere large, O(v—"/?). 
Then, with the usual Prandtl approximation' of neglecting all terms other than those 
of first order magnitude, Eqs. (3.1), (3.2) reduce to 


0g 1 oP 4 ) 0*q a4e 7 (4.1) 
ae aks ie ee ee » Re esd ? 
15s p : On on. 

06 1 0d 

eis (4.2) 

Os p On 


the equation of continuity remaining as in (3.3). 

Further, Eq. (4.2) can generally be neglected, since the total change of pressure 
along a normal must be O(6?). The derivatives of » with respect to s may accordingly 
be replaced by the total differentials. It will also be legitimate, on the basis of (b) 
and (c) above, to regard s as measured along the surface, and m along normals to the 


surface. 
Substitute, now, (3.3) in (4.1) and integrate with respect to m. Then 


. = dp oq 
rf g?—dn+n—+C = (ut pe) —> (4.3) 
o On ds On 


C being the integration constant. 
For brevity, write r for the total shear stress (u-++pe)(0g/0n). Then C is equal to 


the surface value To, and (4.3)-may be expressed as 


T p ” 00 
jatattont+—f q’ — dn, (4.4) 
Tov 0 dn 


TO 


with w;= d6dp/tods, n=n/6. Hence, from (4.4) 


a6 ¢ ) 
— = — [|— — 5}, 
dn ~— pg? \An 


or putting X*=79/pqj, where q is the velocity at the outer limit of the boundary layer, 
06 xX? 0 
~ -—— (z- ws), (4.5) 
on ( q ) dn 
q1 


* It is usual in the dimensional analysis.of the equations of motion to regard 6 as O (v”/?), which is 
the case if the viscous terms are taken to be of the same order as the inertia terms. However, it is known 
in the case of the turbulent boundary layer that 6 is proportional to v5. Cf. S. Goldstein, Modern develop- 
ments in fluid dynamics, vol. 2, University Press, Oxford, 1938, p. 362. 

5 L. Prandtl, Uber Fliussigkeitsbeweggungen bei sehr kleiner Reibung, Verhand. des dritten internat. 
Math.-Kongress, Heidelberg, pp. 484-491 (1904). 
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so that 
ey yo 
6= x f ———— | — — os }dn + const. (4.6) 
A i 
qi 
5. Equation of the boundary layer parameter H. If the constant of integration 
in (4.6) is chosen so that @ is zero at the surface (n=0), it follows that —@ will be 
equal to the angle between the tangent of the streamline at the point P (Fig. 1) and 
the tangent to the surface where the orthogonal through P strikes the surface; or in 
less precise terms, —@ measures the divergence of the mean flow with respect to the 
surface. 


Let 6; be the value of @ at »=1, s; the distance along the outer edge of the bound- 
ary layer, and let the value yj; of the stream function y be defined as 


5 
V1 -f gdn. (5.1) 
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we have, to the order of accuracy of the boundary layer equations (si= s), 
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Now introduce the displacement length 6, which, in terms of 6, may be written non- 


dimensionally as 
s 1 
Me -f (: ‘) dn. (5.3) 
6 0 q1 


With respect to (5.1) we have, therefore, 
¥1 = qi(6 — 41), (5.4) 
so that from (5.2) and (5.4) 
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and, on differentiating (5.5) with respect to s and rearranging, 
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There is also the momentum equation which takes the well known form 


dd v dq To 
— —(H+ 2) = — = X’, (5.7) 
ds qi ds pq? 


where the momentum length @ is defined as 


8 
o-f £(1-*) an, 
0 41 q1 


ds ds\8/ #8\ds_ 8 ds/’ 


dH db dd 
aes 
ds ds ds 


and H = 6,/8. Hence 


and therefore 


(5.8) 


Now substitute for d6,/ds from (5.6) and for dd/ds from (5.7). It is then found that 
Eqs: (5.6), (5.7) may be written alternatively as 


dH 
—- — X*[ws + H{1+ w9(H + 1)}] — A, (S.9) 
$s 
dd 
ea X?[1 + we(H + 2)], (S. 10) 
as 


in which w,= ddp/rods and wy=0dp/Tods. 

6. Discussion of Eqs. (5.9), (5.10) in relation to the empirical formulae of von 
Doenhoff, Tetervin and Garner. It is shown in Part 2, Sec. 12, that Eq. (5.9) can be 
reduced to the form 

dH 
o—— = Fy(w9, Rs, H), (6. 1) 

ds 
R, being the Reynolds number giW/v. This is of interest in that von Doenhoff and 
Tetervin' recently introduced an empirical equation in terms of wy and H, that is to 


say, 


dH 
3 — = F,(we, H), (6.2) 
ds 
though analytically their expression differs appreciably from the theoretical relation. 
Following a similar procedure, Garner? has derived another empirical equation which, 
functionally, is of precisely the same form as the theoretical result (6.1), but which 
otherwise bears a close resemblance. to the formula of von Doenhoff and Tetervin. 
Garner’s solution, therefore, is of particular interest in relation to Eq. (5.9), and it 
may be useful to summarize briefly these empirical developments. We will begin by 
considering the original work of von Doenhoff and Tetervin. 

They observe first that the velocity profile determines H, and point out that the 
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converse statement cannot be proved theoretically. They accordingly proceed to sub- 
ject the hypothesis to experimental test, and, from a considerable volume of experi- 
mental data, show that there is very convincing evidence to support the assumption 
that H uniquely defines the velocity profile. On the strength of this conclusion that 
the distribution is a uni-parametric function, they then consider how the external 
forces acting on the boundary layer are related to H. The argument is advanced that 
the rate of change of H, rather than H itself, is the determining factor, and it is further 
pointed out that this assumption has the desirable effect of connecting conditions 
downstream from a point with those upstream of the point. The problem which von 
Doenhoff and Tetervin investigate, therefore, is the degree of correlation between 
dH/ds, the pressure gradient dp/ds and the surface friction 79. In the first instance 
they attempt to establish a relation between the gradients of H and p when expressed 
non-dimensionally as* #dH/ds and —dd(p/4pq;)/ds, thus leaving the frictional term 
as an independent entity. From their analysis, however, the authors conclude that 
there is no general relationship of this kind, a systematic variation with Reynolds 
number being noted. They then consider the ratio of the dimensionless pressure gradi- 
ent, given above, to the friction intensity which they write in the.non-dimensional 
form To/pqj. This ratio, it will be seen, is equal to —28dp/rods. It is therefore propor- 
tional to ws. In determining the above quantity, von Doenhoff and Tetervin tenta- 
tively assume the flat plate skin friction law, as given by Squire and Young, irrespec- 
tive of the pressure gradient. Allowing for the fact that dH/ds and dp/ds were ob- 
tained by graphical differentiation of the experimental results, they conclude that, at 
a constant value of H, there is an approximately linear relationship between ddH/ds 
and #dp/rods, and consequently arrive at the general result indicated by Eq. (6.2). 
Finally, by analyzing their data at a number of prescribed values of H, they are able 
to formulate an arbitrary, analytical expression for (6.2) which, in terms of the pres- 
ent notation, is given in section 13, Eq. (13.1). 

Garner follows essentially the line of development instigated by von Doenhoff 
and Tetervin. He expresses the momentum equation in the form first used by Ho- 
warth,’ and, as a result, prefers a power law for the skin friction.* He also takes into 
account the effect of pressure gradients on the skin friction, whereas von Doenhoff 
and Tetervin are content in general with the plane flow approximation. Using © and 
I’ as basic variables, where** O=0R;", and T'=Odq:/qids, Garner then obtains an 
equation for @dH/ds which is analytically similar to that of ref. 1. In the existing 
notation it may be expressed, however, in the alternative form given in section 13, 
Eq. (13.2), from which it will be seen, since X is a function of Ry, that it isentirely 
consistent with Eq. (6.1). 


* The presence of the negative sign will be clear when it is recalled that von Doenhoff and Tetervin 
adopt the alternative form #dg/gds, where, in their notation, g is the dynamic pressure at the outer limit 
of the boundary layer. 

** The index is subsequently taken as  =6 in conformity with the choice of Falkner’s equation for 
the skin friction. 

6 H. B. Squire and A. D. Young, The calculation of the profile drag of aerofoils, Tech. Rep. of the 
Aeron. Res. Comm., R & M. No. 1838 (1938). 

7 L. Howarth, The theoretical determination of the lift coefficient for a thin elliptic cylinder, Proc. Roy. 


Soc. (A) 149, 574-575 (1935). 
8 V. M. Falkner, A new law for calculating drag, Airc. Eng., 15, 65-69 (1943). 
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The empirical equation for ’dH/ds, the momentum equation and a relation for the 
skin friction are then sufficient to enable the development of the turbulent boundary 
layer to be analyzed. The solution depends on a step-by-step integration which Garner 
elaborates in great detail by the finite difference calculus. It assumes that th€ static 
pressure distribution is given, and that the initial values of 8 and H are known, e.g. 
at transition. This problem is discussed thoroughly in both refs. 1 and 2, and need not 
concern us further. The integration then yields 3 and H at the beginning of each in- 
terval. Hence, we seek a form of Eq. (5.9) purely in terms of 3, H and p. This is the 
subject of Part 2. 


Part 2. 


7. Approximations regarding the turbulent and laminar layers. Boundary condi- 
tions. It is inevitable that a theoretical approach to the-problem of the turbulent 
boundary layer must be appreciably more complex than the simple, empirical treat- 
ment of von Doenhoff and Tetervin. The difficulty of solving the fundamental equa- 
tions of motion, even when the flow is laminar, is here increased by the fact that we 
have to consider an “apparent viscosity” which varies from point to point in the fluid. 
In order to deal with this feature of the flow, some kind of turbulent mechanism must 
be specified and incorporated into the basic equations. The very considerable diffi- 
culties of a rigorous treatment of this aspect of the problem are well known, and even 
an approximate and much simplified theory, as adopted in the present paper, leads 
to a rather complicated solution which does not yield an analytical function for 
idH/ds. Indeed, it has not been possible to present the ‘arguments advanced in the 
following pages in the form of a unified theory. Rather the investigation has been 
divided into a series of associated problems which are considered individually, and it 
is then shown how the results may be combined in a numerical solution to provide 
the data for the step-by-step integration referred to in section 6. As a consequence, 
it will probably assist the reader to give a statement as to the procedure now to be 
followed, and the nature of the approximations involved. 

In the first place it will be assumed that the flow may be sub-divided into a tur- 
bulent layer, where the effects of fluid viscosity are negligible, and a thin surface 
layer (laminar sub-layer) where viscosity is predominant, the transition from thé one 
flow to the other being regarded as occurring instantaneously, i.e. in zero length 7. 
Neglect of vy when the motion is fully turbulent implies that €>v (see Eq. 4.3), an as- 
sumption which is valid for sufficiently large Reynolds numbers, and one commonly 
made in considering this region.® It will also be convenient at first to take ws as in- 
dependent variable, and we will subsequently derive a method of transforming from 
ws to wy which is the variable required in performing the numerical integration of 
Eqs. (5.9), (5.10). 

On the above basis, we shall then proceed to our first objective, namely the estab- 
lishment of the velocity distribution in the turbulent and sub-laminar layers, the two 
cases being considered independently, but in such a manner that the essential bound- 
ary conditions are satisfied, and further, that the velocity is continuous at the transi- 
tion from the laminar to turbulent states. Secondly, the important quantity 6; of 
Eq. (5.9) will be investigated. This, it will be shown, depends primarily on the condi- 





® See, e.g., S. Goldstein (editor), Modern developments in fluid dynamics, vol. 2, University Press, 
Oxford, 1938, pp. 331-332. 
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tions in the laminar sub-layer, which cannot be entirely neglected for this reason. 
To a less, but appreciable degree, it also affects the magnitude of the parameter H. 
Finally, we develop a method whereby Eq. (5.9) may be evaluated numerically in 
terms*of ?, H and p (or qi, which is related to p by Bernoulli’s theorem). 

Turning, now, to the approximate nature of the analysis, we have first the simple 
theory of diffusion upon which the flow in the turbulent layer is based. It is assumed 
that the intensity components of the turbulent fluctuations do not differ appreciably 
(or at any rate are proportional to one another), and that the scale may be sufficiently 
represented by a mean length which is a function of the space co-ordinates. Although 
the subsequent development of this simplified conception of the turbulent mechanism 
near a surface in no way depends on a physical model, it is mathematically equivalent, 
nevertheless, to Prandtl’s momentum transfer theory,'® and is further closely con- 
nected with von Kaérman’s similarity theory"' which appears as a special case in the 
present treatment. 

Without entering into the details of the argument, it is finally shown that the tur- 
bulent velocity distribution may be expressed in terms of two functions f and g, and 


the surface friction X, namely, 
) ( q ) X 
dn \ qu ciel foul 


” 
Afro? — -{ gf !*dn a const., 
0 





where 


X being a function of the correlation between the longitudinal and lateral velocity 
fluctuations, and the scale of the turbulence. 

Strictly, when the natural viscosity of the fluid can not be ignored, f and g are de- 
pendent on w;, 7 and the Reynolds number, but, when the flow is fully trubulent, we 
neglect, in accordance with our initial assumptions, the viscous terms, and treat both 
fand g purely as functions of w; and yn. The stress functionf may then be equated to fr, 
where fr denotes the component due to the Reynolds shear stress. Hence, assuming 
for the moment that X is known, the problem of calculating the velocity distribu- 
tion in the turbulent layer reduces to the determination of fz and g. To this end, we 
first make use of an argument advanced in section 8 (a), namely that approximate 
solutions for fg and g are sufficient provided Afg/? is itself correctly established. 
Accordingly, we express both fz and g as a power series in 7 with coefficients which are 
to be purely functions of ws. In the case of fe the coefficients are chosen to satisfy all 
the known boundary conditions, with the exception of terms dependent on viscosity 
which are regarded as referring to the laminar sub-layer only. This leads to a solution 
in quartic form. In addition to the series for g, two particular solutions for \ are con- 
sidered, viz. (a) when g=function of ws only = —A, (b) when \=function of ws only 
=),. It is then suggested that (a) holds near the surface, while (b) is applicable to- 


10 L. Prandtl, Bericht tiber Untersuchungen zur ausgebildeten Turbulenz, Zeitschr. angew. Math. Mech. 


5, 136-139 (1925). 
11 T, y. Kaérmén, Mechanische Ahnlichkeit und Turbulenz, Nachr. Ges. Wiss. Géttingen, Math-Phys. 


Klasse, 58-76 (1930). 
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wards the outer edge of the boundary layer. By identifying A: with A when n=1, 
and assuming that 0A/0n is then zero, a solution for A, in terms of A and \y may 
then be obtained, if g is expressed as a quadratic, to satisfy condition (a) when 7 =0, 
and condition (b) when n=1. The resulting distribution of \ for intermediate values 
of 7 is then regarded provisionally as a valid approximation. Further, when w;=0, it 
appears that the corresponding value of A (designated Ao) is then identical to. von 
K4rm4n’s constant K of the similarity theory."! Its value is therefore known. In 
addition, we obtain immediately an integral relation for A/Ao purely in terms of fr. 
Hence, apart from X, A; is the only remaining unknown. As a first approximation, 
but nevertheless one which appears to be well substantiated by experiment, we as- 
sume that A, is independent of ws. Like Ao, it must then be regarded as a fundamental 
constant to be determined experimentally, the value in the present instance being 
calculated to give the best agreement between the theoretical and experimental veloc- 
ity profiles for flow in parallel wall channels. Finally, neglecting second order terms, 
we obtain a very simple relation for the variation of X with w;, namely X/Xo=A/Ao, 
where Xo refers to the condition w;=0, and is therefore calculable from either von 
K4rmA4n’s logarithmic skin friction equation,"-'* or from a power law such as Falkner 
has published.* We thus establish a general solution for the velocity distribution in 
the turbulent part of the boundary layer. 

As regards the treatment of the laminar sub-layer, little need be said. It is essen- 
tially in the nature of a linear (double link) interpolation which satisfies the main 
wall condition, and preserves continuity with the turbulent flow solution at the point 
of transition. It makes no reference, therefore, to the equations of motion. On the 
other hand, the approach is justified on the grounds that, when the Reynolds number 
is moderate or large, the laminar layer is quite thin, and further, that the velocity 
distribution is then mainly linear throughout the region concerned. 

The solution for 6, however, requires a little more attention. As already pointed 
out, it is primarily dependent on conditions in the laminar layer; it is also largely 
determined by the stress function. Hence, the relation for fg considered in the study 
of the turbulent velocity profile is certainly no longer tenable near the surface. We 
accordingly develop the stress function, which must now be written as f, as a new 
power series, applicable for small values of 7 only, and including the effects of viscos- 
ity. 

The procedure for determining the coefficients of the series is similar to that 
adopted in dealing with the turbulent layer, except that in this case the boundary 
conditions are restricted to the surface. The series is also limited to a quartic in view 
of the fact that n is to be regarded as quite small; at the same time it is desirable on 
account of the increasingly involved character of the coefficients relating to higher 
powers of 7. The distribution of f in the outer part of the boundary layer is then 
represented by a second, independent series which satisfies the conditions at »=1, 
and is continuous with the distribution of f near the wall. The fact that this “inner and 
outer” solution is not altogether consistent with the previous solution for f, when 
viscosity was entirely neglected, is regarded as unimportant, for in the first place both 
approaches are only approximate, and secondly, so far as 6; is concerned, it is the dis- 
tribution of f in the laminar sub-layer which is critical, whilst the relatively large 





12 T. vy. Kérman, Turbulence and skin friction, J. Aer. Sci. 1, 1-20 (1934). 
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error which probably arises in the turbulent region has only slight significance. In the 
former case there is reason to believe that, when viscosity is included, the treatment 
leads to results of good accuracy. 

Finally, having established @;, numerical estimates of ’dH/ds may be made from 
Eq. (5.9) in terms of w;. An additional equation relating ws; and wy is also given. Hence 
Eq. (5.9) may be integrated step-by-step according to ref. (1) or (2). 

8. (a) Turbulent diffusion in the boundary layer. The first step towards the trans- 
formation of Eq. (5.9) is to develop a theory of turbulent diffusion, valid for the 
boundary layer. This is greatly complicated by the fact that the flow is anisotropic. 
Nevertheless, an approximate treatment, which has yielded valuable results in the 
case of turbulent flow over plane surfaces and in pipes, rests on the assumption that 
the intensity components 7/4”, ~/v” do not differ appreciably, and further that the 
correlation may still be represented by a single scalar function. Extending this hy- 
pothesis to the case of flow in the presence of pressure gradients generally, we write 


e— l/7?, (8. 1) 


where / is the length defining the average scale of the turbulence at any point. To the 
order of accuracy of the boundary layer equation pedg/dn is equal to the original 
Reynolds stress —pu’v’. Denote this component by rr. Then 
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where the correlation coefficient Ry» is here defined as 
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since, by hypothesis, y/t=y’?, 
From (8.2) we then have 
= € 0g 
vy? =——_ — (8.3) 
Rae On 
and combining (8.3) with (8.1), 
2 dq 
€¢€ = —= 9 (8.4) 
Ri» On 


so that the Reynolds stress may be written non-dimensionally as 
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the function \=R;,*(1/4) being, in general, a function of 7, ws and the Reynolds num- 
ber, since the turbulent mechanism is influenced by viscosity near the surface. It fol- 


lows that 
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We now seek to express \ in terms of the stress function fr. First differentiate (8.6) 
with respect to 7. Then 


Or 0° 0 1 ¢@ 
a [= of a =o. (8.7) 
On On? On 2fr an 
Without loss of generality we may also write 
re) 3? 
beget I aes (8.8) 
On On? 


g being also a function ef 7, ws and the Reynolds number. Hence, combining (8.7), 
(8.8) 


————=-y, (8.9) 


which on integration yields 


, 7 _ 
be~e | f gfe dn + const. ' (8. 10) 
0 


- 
or writing J; = — {¢gfx/*dn+const., 
d= fe Ih. (8. 11) 


Ignoring for the present the integration constant in (8.10), consider, now, the solution 
when g is a function only of w; and Reynolds number, so that with respect to 7 


g = constant = — A. (8.12) 
Then 
” — 
A= Afe f fr’ dn, (8. 13) 
0 
for which there is a maximum at 7=m such that 
_2 { Ofr “ 1 
fren (=) f fr dn = —2. (8. 14) 
on hl 0 
Again, the solution which makes 
A = constant = Ay (8.15) 


for all values of 7 is given by 


Pye Kee (8. 16) 


Now assign to \; the value of \ when 7=1. If, further, we make the reasonable as- 
sumption 0A/8n-0 as n—1, then a solution for \ is obtained, which satisfies the above 
conditions, by equating (8.11) to 4: when 7=m, with A=A, in the range m<7<1. 
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Whether or not this solution for \ leads to boundary layer characteristics in accord- 
ance with fact can only be decided by test, but for the present we will regard it provi- 
sionally as satisfactory. Hence the necessary and sufficient conditions are 


n = 0, A = 0, g=-A | 
1 ofr r (8.17) 

7 = 1, A=A,, Cee (=) A1| 

2frm \On/ mn) 


In general write 


g=- a(1 +> Ar") : (8. 18) 


n=1 
where A, are arbitrary constants. Then (8.17) will be satisfied by the series 
g = — A(1 + Ain + Azn’), (8.19) 


so that from (8.10) 
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Go = fim f Sr dn, G;= fie, { nfr dn, 
0 0 
: 1/2 "=i 
Ge = Sem n°fr dn, 
0 
and 


Ii 


=" rs] m™ _wye 
Ho =1+ 4am. (=) f fr dn, 
On 1 0 


-12 [ ofr “ —3 
Hy = + 3frn (=) f fr dn, 
"1 0 


On 


2 -12 [ ofr " —1/2 
He = 71 + frm (=) f n’fr “an. 
On 1 0 


<) 


Then (8.14) makes 
Ho = 0, (8.22) 


and, to satisfy (8.17), we have from (8.20), (8.21) 
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A 
Ap (8.23) 


AyAy + Hohe = 0 


Ai 
He (— _ Gv) 
A 


Go + GiA1 + Gedo 


Therefore, 





Ai = , (8.24) 
(Giz — G2H;) 
Ay, 
Ao =_— z (8.25) 


Hence, by the definition of h, when 0<y<m, 
’ 2 
i= af (1 + Ayn + Aen) fr” dn + const.; (8.26) 
0 


and when 7,<7<1, 
= ifr + const. (8.27) 


It follows from Eqs. (8.6), (8.11) that 


= () _ (8.28) 
an\as Th 


whence, integrating along normals from the outer edge of the boundary layer towards 
the surface, and allowing for the condition g—q, as n—6, 


£ =1+/2X, (8.29) 
q1 f 


where I,=/?I7'dn. Hence, the velocity profile is determined primarily, not by the 
individual values of fr and X, but only by the product Afz"/?, for which there will be 
a unique distribution across any section of the boundary layer, depending on the flow 
conditions. We may therefore seek approximate solutions, of an arbitrary character, 
for fg and X, provided that, when combined in the above manner, they yield a solution 
of (8.29) in accordance with the physical facts. 

In the next section we shall make use of this argument to obtain a general solution 
of the velocity distribution in the turbulent layer. 

8. (b) Velocity distribution in the turbulent layer. It has been shown experimen- 
tally'*.4- that, for moderate or large Reynolds numbers, fully developed turbulent 
flow in the absence of a pressure gradient, or when the gradient is small (as in pipe 


13H. Darcy, Recherches expérimentales relatives au mouvement de l'eau dans les tuyaux, Mém. Sav. 


Etrang. 15, 141-403 (1858). 
“4 T. Stanton, The mechanical viscosity of fluids, Proc. Roy. Soc. (A) 85, 366-376, (1911). 
6 W. Fritsch, Der Einfluss der Wandrauhigkeitsverteilung in Rinnen, Zeitschr. angew. Math. Mech. 8, 


199-216 (1928). 





196 W. S. COLEMAN [Vol. V, No. 2 


flow), does not depend on viscosity. This fact is expressed quantitatively by the well 
known velocity defect law which, in the present notation, becomes* 


1 


x 


(8.30) 


te oes 


“) = ¢(n), 


71 
or, extending (8.30) to include the case when the pressure gradient is not negligibly 
small (w;+0), 


— Ts 


1 q 
~hea(t= *) = $(ws, 1). (8.31) 
A q1 
Hence, 
a, 1 fe ase 
ee 2. 8.32 
; 7, . d (ws, 7 ( ) 


But, from the preceding section, we have obtained X, and therefore J,, in terms of fr 
and the variables A and \; which are purely functions of ws. Hence, it follows that in 
the turbulent region, for which Eqs. (8.31), (8.32) are valid, fp must be a function of 
ws; and 7 only. Moreover, we may then write f—fr, and therefore f is also simply a 
function of w; and 7. 

We will now obtain an expression for J;, which satisfies the turbulent velocity 
distribution generally, by the argument advanced in the last section. This, it will be 
recalled, depends only on the approximate determination of both fr and X, provided 
that Afz"/? is correctly established. First, therefore, express fr in an arbitrary form. 
The problem then reduces to the determination of A and \; so as to satisfy the above 
conditions for J;. Accordingly, proceeding to the approximate development of fr, 
we note first that Eq. (4.4) has the alternative form 

(+) dn. 
q1 


Ae 


Successive differentiation of (8.33) also gives 


re] 
Os 





5 
X? 


f=1+0m+ (8.33) 








of 5 (£) ce) (=) 

a oe wee b=, eee 2, 

On X* \qs os\q 

ae ori riciarereses 

On? ‘ XxX? On q1 Os q1 q1 Os On q1 , 

a3 67a /q\a/q a/q\a a a 3 

4 SEM rtQEL+(92 (9) 
On? De On? q1 Os q1 On 71 Os On qi q1 Os On? q1 
etc. 





Consider, next, the known boundary conditions. They are: 


* This equation is commonly written in the form (U.—U)/Ur=f(y/h). 








nN 
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qa rs] rs] 6 
it ee = (+) a ee (4) = 2° y2 = R,x?, 
91 Os\Q On \ 1 v 
d i ¢ ri) 
oi, Dag - gare —(+)=0 
q1 ds p ds . On \ qi 
Hence, for f and the known derivatives, we have 
fe) 3? 
n = 0, f=1, ll by het iP 
On On? 
38 dX | 
os = axRis | + — @) ‘ 
On?® d. qi ds 
and when 
of 


—= 0. 
on 

These conditions must be satisfied in addition to that discussed earlier, namely, for 
fully developed turbulence, the flow is independent of the Reynolds number when 
moderate or large. It also follows that fp may then be taken equal to f. Express, 
therefore, fg as a power series in 7, with coefficients which are purely functions of a. 
If, in addition, fr is to satisfy the boundary conditions, we must ignore those in which 


the Reynolds number occurs explicitly, the argument implying that such terms are 
only important in the laminar sub-layer. This leaves five boundary conditions to be 


satisfied. Accordingly, we write* 
‘fr = Ao + Arn + Aan? + Agn® + Aan’. (8.34) 
Hence, 
Ao = 1, A, = ws, A; = 0, (8.35) 
and 
Ao +Ai1+Ae+Azs+4,=0 
A, + 2A2+ 3A3 + 444 = of 
which on solution yield 
A; 
A, 


ay | (8.36) 


3 + 2w; 


We must now consider the boundary functions A and \,. First, with respect to Ai, 
let us provisionally assume that it has a unique value for the fully developed turbu- 
lent boundary layer (i.e. it is constant for all positive values of w;). It then follows 


from (8.13) that 


* The above series, to fit the same boundary conditions, was first suggested by Fediaevsky. See 
K. Fediaevsky, Turbulent boundary layer of an aerofoil, J. Aer. Sci. 4, 491-498 (1937). 
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(8.37) 


the suffix 0 to A and the bracket of the numerator on the right hand side of the equa- 
tion denoting the condition w;=0. Hence, it remains to determine A» and \y. From 


Eqs. (8.8), (8.19) we see that, when 7 =0, 
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which is recognized as aform of von Kaérm4n’s equation for the length/ in the similarity 
theory." Accordingly, when w;=0, A must be equal to the Kérm4n constant K whose 
value, from the logarithmic skin friction law for plane flow (see section 8 (c), Eq. 8.47) 
is 0.392. Thus, when 


ws; = 0, A = Ao = K = 0.392. 


Finally, if \; is to be independent of ws, we may determine its value for any arbitrary 
condition of turbulent flow, e.g. for the case of plane flow. Like Ao, it must also be 
regarded as a constant, only ascertainable numerically by experiment. We will there- 
fore assign to \, that value which gives the best theoretical agreement with the data'®-!7 
of Fig. 2. This leads to \; = 0.14, when it will be seen that the distribution of I: (ws =0) 
calculated from the numerical integration of Eqs. (8.26), (8.27), in conjunction with 
Eqs. (8.34), (8.35), (8.36), agrees very well with the experimental curves for all values 
of » outside the laminar sub-layer. 

The validity of the assumption A; =constant, and hence of Eq. (8.37), is discussed 
in section 8 (c). Further evidence appears from the general comparison of the theoreti- 
cal and experimental velocity profiles in section 10. 

8. (c) Consideration of the general skin friction law (w; 20). The present develop- 
ment of the theory of turbulent flow near a surface depends essentially on the three 
boundary terms A, A; and X, of which A and ),; have already been considered. As 
regards X, which determines the surface value of the shear stress, or skin friction 
intensity, previous investigators have mainly been content to ignore the effect of 
pressure gradients on X, and to take one of the well attested laws strictly applicable 
for plane flow only. We shall now show, however, that a simple, approximate relation 
may be derived to account for the variation of X with ,, i.e. for the influence of pres- 
sure gradients on skin friction. This leads to the consideration of the flow conditions 
very near the surface, namely for very small values of 7. Eq. (8.34) then reduces to 


fr =1+ wm, (8.39) 
and, equally, Eq. (8.19) approximates to 
g=-—A. (8.40) 
Substituting, now, (8.39), (8.40) in (8.26) and integrating, we find 
2A 





(V1 + wan — 1), (8.41) 


2 


I, = 
Ws 
the constant of integration ‘satisfying the condition J,=7=0. Or, since wn is small 
in relation to unity, (8.41) may be expanded to give simply 


I; = An. (8.42) 


Hence, by definition, 


16 F. Dénsch, Divergente und konvergente turbulente Strimung mit kleinen Offnungswinkeln, For- 


schungsber. Ver. Deutsch. Ing., p. 282 (1926). 
17 J. Nikuradse, Untersuchungen iber die Strémungen des Wassers in konvergenten und divergenten 


Kanilen, Forschungsber. Ver. Deutsch. Ing., p. 289 (1929). 
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1 1 1 
I, = — dy = — f — dn, 
q; A | 


1 
— log. » + const., 
A 


ho + hy logio n, (8.43) 


where ho=constant for any given value of ws, and 4; =2.3026/A. Owing to the fact 
that viscosity has been neglected, Eq. (8.43) leads to the well known result that 
g=-—« when 7=0, or, alternatively, that g=0 at some small distance mo from the 
surface, mo being clearly of the order of thickness of the laminar sub-layer. As shown in 
the appendix, the original boundary layer equation (4.3) indicates that mo is propor- 
portional to 1/R;X, and we therefore write 

Q 


= 8.44 
R;X ( 





? 


70 


where 2, the factor of proportionality, is clearly a function of ws. It is also immediately 
apparent from (8.29) that, for the approximate condition n=, g=0, 


1 
(2) = — +: (8.45) 
Hence, combining Eqs. (8.43), (8.44), (8.45), 
1 
XY == (ho + hy log 2) a hy log R;X, (8.46) 


which has the form of von Kaérm4n’s skin friction law for plane flow.'':'* Thus, when 
w;=0, Eq. (8.46) must be consistent with von Ka4érm4n’s semi-empirical relation 


1 

— = 3.60 + 4.15 logy Ric; , (8.47) 
cyl? 

where cy/?=+/2X and, by identifying (8.46) with (8.47), it follows, when fy and fh; are 
determined from Egg. (8.34), (8.35), (8.36) for w;=0, that 2 has the value 0.131 which 
compares with values given by Prandtl'* varying between 0.089 and 0.111. Eq. (8.46) 


then reduces to 


ee + 5.88 log RsXo, (8.48) 
0 

and it follows immediately from the relation K =Ap =2.3026/(M)o that Ao=0.392. 

For the case when w; 0, it is necessary to know Q as a function of w, in order to 
obtain X from Eq. (8.46). The difficulty of establishing this relation theoretically may 
be avoided, however, in the following approximate treatment, which appears to be 
reasonably valid up to values of w; at which separation is imminent. 

On splitting the product of Reynolds number and skin friction coefficient in (8.46), 
we have the alternative arrangement 


18 L. Prandtl in W. F. Durand, Aerodynamic theory, vol. 3, J. Springer Berlin, p. 140 (1935). 
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1 
> = hy, log Rs + hy log X — (ho + hy log Q), (8.49) 


or by writing 
hy log X — (ho + hy log Q) = ohy, 


where @ is also a function of ws, 
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X = h,(log Rs + a). (8.50) 
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Fic. 3. Comparison of X/Xo and A/Ao. 


Consider, now, the variation of X with w; when R; remains constant. Then 
X h{ (log Rs + o’) 


Bw ; (8.51) 
I h,(log R; + o) 





where i;, X, ¢ correspond with ws, and hy’, X’, o’ with w,’. 

It further appears from experiment that @ varies only slightly, and is small in 
relation to log Rs, provided R; is moderate or large. Hence, for the ratio (8.51) no 
serious error is introduced in neglecting o, and we then have the simple relation 


eee lie 
pa sents ta oee'g (8.52) 


xX’ hy A’ 





which, when w,’=0, becomes 


x 
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A/Ao being given by (8.37) and Xo by (8.48). 

Evidence!®”° as to the validity of (8.53) is shown in Fig. 3. 

From the very limited data at present available, it would appear that (8.53) is a 
reasonably good approximation in the range 0 <w;<50, but clearly a good deal more 
experimental information is required before the precise significance of Eq. (8.53) can 
be ascertained. 

9. Laminar sub-layer. The solution of Eqs. (5.9), (5.10) depends on a knowledge 
of the momentum thickness, and the parameters H and 4, each of which is a function 
of the velocity profile. If, however, the laminar sub-layer is neglected in evaluating 
these quantities, appreciable errors will arise, and it is necessary, therefore, to consider 
the flow in this region as well as in the fully developed turbulent part of the boundary 
layer. Before doing so, it will be convenient to express the displacement and momentum 
lengths in the non-dimensional forms 


A 6 1 
—=f (1-2) dy = ax; (9.1) 
6 0 q1 


<= 








! é 'q q 
\ —= f —{1- £) dn = aX — BX?, (9.2) 
\ 6 i] q1 71 
ictal where, from Eq. (8.29), a= —J$I2dn, B= fp I3dn. Con- 
‘\ sequently, 
\ 
d a? (/H — 1 
ons Sih ee ) (9.3) 
, 6 B\ #? 
Te] 1 
6 ie (9.4) 
-I, 8 
1——%X 


a 


Turning our attention, now, to the lami- 
nar sub-layer, we note first that, for mod- 
erate or large Reynolds numbers, it is 
quite thin. Since, also, the velocity dis- 
tribution in this region is mainly linear, 
a good approximation to the correspond- 
ing distribution of Jz may be represented 
by the discontinuous curve ABC (Fig. 4). 
The characteristics of his curve will be 
Fic. 4. Distribution of J; in laminar sub-layer. determined by satisfying the essential 











19 A. Fage and V. M. Falkner, An experimental determination of the intensity of friction on the surface 
of an aerofoil, Proc. Roy. Soc. (A) 129, 378-410 (1930) ; also Tech. Rep. of the Aeron. Res. Comm., R & M. 
No. 1315, 1, 117-140 (1930). 

20 G. B. Schubauer, Airflow in the boundary layer of an elliptic cylinder, Twenty-fifth Ann. Rep. Nat. 
Ad. Comm. for Aeron. Rep. No. 652, 207-226 (1939). 
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wall conditions, and by preserving continuity with the solution (v neglected) where 
the laminar flow merges with the turbulent flow. For 


0O<n<7, 
let 
Ts = jo + jim (9.5) 
and for 
n<n<(n' +7"), 
Tz = jo + (j1 —J2)n’ + jan, (9.6) 


where jo, ji, je are coefficients satisfying the necessary boundary conditions, and 9’, n’’ 
are defined in Fig. 4. Then, if ; is the nominal thickness of the laminar sub-layer, 


vm 


AR;X 


(9.7) 





m= +7" = 


according to the appendix. 
Also, from Eqs. (8.45), (9.5), we have, when the approximate condition 7 =p is 
replaced by the correct condition 7=0, 


1 
jo = , a (9.8) 


and, from the differentiation of (9.5) 


Gat 
J = — = —— e 
On / <0 11] y=0 
But, according to Eq. (8.28) 


Cr). x la 


which, from the boundary condition 


gives 


1f9/4q 
nr=— E £) = RX. (9.9) 
X Ldn \ qi/ Jy=0 


Finally, if the distribution of J; in the laminar layer is to be continuous with that in 
the turbulent region, we have from (9.6) 


2,0), om 
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where (J1)»,, (2), are the values of J; and J, at 7 = given by the solution of sections 
8 (a), 8 (b) when » is neglected. 
It follows from (9.7), (9.9) that 


/. 
Vm 





m=y +9” = — 
An 
and hence, from the relation 


jo + jin’ = (T2)n, ~ jan”, 





that 
(oes pa a (9. 12) 
(j1 — je) 
Hence, by definition, 
n’ m1 1 
a= — f° Got jmin— f" bint Ga onl + jon}dn — fi Tad 
0 7’ 1 
1 
= = jon a 3(jin”? +} 2j1n'n” + jon”) -f Iedn, (9.13) 
"1 
and 
7! ml 1 2 
B -f (jo + jrn)*dn +f {jo + (j1 — j2)n’ + jan} %dn +f I edn, 
0 9’ "1 
= jon, + jrljo + jin’’)n’? + jojo + jin’)n’” 
1 
2 a, .2 / 2 ai . 
+ 2jojrn'n” + 3(jin’* + jon’”*) +f I dn. (9.14) 
v1 


10. Velocity profiles. Comparison of theory and experiment. Von Doénhoff and 
Tetervin' have analyzed a large amount of data from velocity measurements in the 
turbulent boundary layer, and it is of interest, therefore, to compare their empirical 
curves with the theoretical conclusions of the present paper. For this purpose it has 
been most convenient to take the actual velocity distribution in the form g/q =f(n/#) 
as a basis of comparison, the experimental data being obtained from cross-plots of 
Fig. 9, of the paper quoted in Footnote 1, where g/q: is plotted with respect to H 
for a series of values of n/3 =constant. 

Fig. 5 gives typical results at w;=30 and for log R;=3.5 and log R;=6.0. In 
making the theoretical calculations, J; for the turbulent layer (v neglected) is first 
determined from the numerical integration of the reciprocals of Eqs. (8.26), (8.27) 
in conjunction with Eqs. (8.34), (8.35), (8.36) and (8.37); X then follows from (8.37), 
(8.48), (8.53), and hence the velocity distribution in the turbulent part of the bound- 
ary layer from (8.29). The integrals a, 8 are next evaluated in the turbulent region 
by numerical integration of Jz and J} respectively, and thence for the laminar layer 
by direct calculation from Eqs. (9.11), (9.12), (9.13), (9.14). Finally, H follows from 
Eq. (9.4), 8/6 from (9.3) and hence n/? from 
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| 3 aap 


The results of Fig. 5, which are representative of a number of such calculations, indi- 
cate satisfactory agreement between theory and experiment up to values of H in the 
region of separation. 

11. Approximate treatment of Eq. (4.6) in the solution of Eq. (5.9). We are now 
in a position to consider the solution of Eq. (5.9). This equation contains the impor- 
tant term 6, which, as will be seen from Eq. (4.6), depends on the ratio of (0f/0n —ws) 
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Fic. 5. Theoretical and experimental velocity profiles. 


to (q/q:)?. Near the surface, within the laminar sub-layer, 06/dn becomes large, and 
it appears that the value of 4; is critically dependent on the distribution of 00/0n in 
this region. Consequently, viscosity is by no means a negligible factor, and it must be 
taken into account, not only with regard to the local velocity distribution, already 
considered in section (9), but also in so far as it affects the stress function f. In this 
respect, the solution is in marked contrast to that of the turbulent, velocity distribu- 
tion, for which f may be equated to fr. Hence, Eq. (8.34) is no longer tenable, and we 
must consider a more accurate form, at least for small values of 7, when the solution 
will depend predominantly on the viscous terms previously ignored? In dealing first 
with the term 6,, we shall accordingly divide the analysis into two parts; (a) thatin 
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which the effects of viscosity are all important, (b) that region outside the laminar 
sub-layer where, again, viscosity may be neglected. As before, let us represent the 
stress function near the surface by a power series, namely, 


f = Bot By t+ Ba? +:-:-, (11.1) 


and solve for the B coefficients from the boundary conditions at »=0. Then, as before 
(see Sec. 8 (b)) 


Bo = Ao = :: B, = A = Ws, B, = Ao = 0, (11.2) 
08 2 dX x d 
B= = (54) J — ee =), 
6 \ dn? 3 L ds ins ds 


2 »fdX dRp X da 
xRe | — I, (11.3) 


oo dR» ds qi ds 


where Ry =qiW/v. 
To a first approximation we now take the relation between X and Ry, as that for a 


flat plate, and determine dX /dR, from Falkner’s power law form.® This in terms of Ry 


may be written 


X = 0.0808R;, (11.4) 

so that 

dX i 

—* Seen”. (11.5) 

dR» 
Also, 

dR ds 8 d 

ar ae - 4/2 +— 4). (11.6) 

ds vy tds qi ds 


Hence, substituting (11.5), (11.6) in (11.3) 


3 P —13/12_3 dd do dqi 6 dqi 2 
— B;= — 0.00673XRz Rs} — + — —I{+ 
2 qi ds 


or, from Eq. (5.10) and the Bernoulli relation 











1 dqx 1 ? 
q1 ds pqi ds” 
3 ai a d 
— Bs = — 0.00673XRy R; | x*{1 ee a 
2 pgix? ds 
6 od 
i ? X‘R;, 
pgix? ds 


so that ° 
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{1 + wo(H + 1)} + wX], 
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2 o —13/12 -y12 
= — 2X3R; | o.0067s ( -) Rs {1+ 0(H + 1)} + wX] . (11.7) 


We also obtain a relation for B, as follows. Equation (4.3) may be written 


tenth ita tee 
xX? 0 q1 on a X?*R; On q1 
+ /2()/2(9). 
X? | On \qi/ | On \ Qu 
which on differentiation with respect to 7 gives 
9/49\|9 (4 
as) lam Geo 


g \? 00 1 8 /q Or 
()2-23@)+> 
qi/ On Rg On? 
re) 3? 
()5(2)-= 
On q1 On? q1 


q1 On 
and when »=0, g=A=0; therefore, from (11.9), 


aaten 
—|—[— = wr", 
R; On? qi n=0 
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On? q1 n=0 














+ 2x? 





or 


We also have from Sec. 8 (b) 


otf 6 [= (2)<(4) +35 (4)s = (4) 
dn* - x? dn? \qis Os \ qi , dn? \qi/ Os On \ qi 
0 0 3 do 8 /@q 
2922) 50) 
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and when 7=0 


(S) Gala ad aes ee 
ant} e-0 X?% Lan?\qi/ as an\gqi)  an\qi/ as an*\qi/ 4d 
From (11.1) there is also the relation 


1 / do 
By 0 (= 
24 \ dn*/ .=0 


Then, from (11.10), (11.11), the first term of (11.12) yields 


(11.8) 


(11.9) 


(11.10) 


(11.11) 


(11.12) 
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1 6 @ da 1 dX X d 
aces oye a8 (+) (+) a — Xi |= reed *) 
8 X* dn? \qi/ os On\ qu qi ds 
3 . 
> Byws, (11.13) 


according to Eq. (11.3). 
Again, for the second term of (11.12), we note first that (11.10) may be written 
alternatively as 
be) Rx dp 
On? 71 7=0 ? ds 


pq 
and therefore 





i = (£)< (4) a. 5? d*p 
8 X? dn\qiJ as an*\qi) 8 — pg? ds? 


or, since w; = ddp/rods = (6/X*pq;)(dp/ds), we have, substituting for 6 in terms of Ws, 


2 € . 1 3: d? db\?2 
LER Esha ase 
8 X? On q1 Os On? q1 8 





ry ry 
Further, 

dp dq 

oe a 

d*p 








dq 2 ] 
= = Soy! + ’ 
ds? “ ( ds ) ale 


and if the Bernoulli velocity distribution is assumed to be linear, which in general is 


a reasonable approximation, 
dp (2) 1 
ds* ds a 


pqi 
Hence, 


1 6 re) q re) 0? q 1 = ee 
—- ee Sage EF ae see oe -erw x. Ryw 3X. 
8 X? dn\qi/ as On? 8 

But from (11.7) 


2 ‘ 3B; 
X*R; = 


d —13/12 1/12 : 
2 | 2.00673 (=) Rs {1+(H+1)+ wX} | 


For brevity write now 
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I 


d —13/12 ~1/12 
= 2 0.00673 (5 R; {1 +. we(H + 1) + wX } ° 


3B, = — oX*R;. (11.14) 


. 


b 


il 


Then the second term of (11.12) reduces to 


1 58 dfq\a &#/q ae 
2, eae ()— itd BS ER ale 3 (11.15) 
8 X? dn\qis das dn? \ qi 8 a 


and adding (11.13), (11.15), 





1 wsX 
By, = 8 bos (1 +- ) . (11.16) 


a 


Higher coefficients become progressively more complex, but in view of the fact that 
the values of 7 concerned are very small, we may, to a sufficient approximation, re- 
strict the series (11.1) to a quartic, in which case we have 


y 
ot is, = 3Bar? + 4Bar' 


dn 
1 wsX 
= nt [1 + - w( 14+ °)y]. (11.17) 
2 a 


This equation, substituted in (4.6), then 
enables the distribution of 0 near the : = 
surface to be calculated. \ 

For the corresponding distribution : 
in the turbulent part of the boundary \ 
layer, we should, to be consistent, cor- \ 
relate Eq. (11.1) with Eq. (8.34), but, \ 
as already pointed out, the value of 6; \ 
depends primarily on the conditions . 
near the surface, and the precise form \ 
of (Of/8n) —w; at larger values of 7 ap- \ 
pears to be less important, provided the \ 





essential conditions relating to f in the 
turbulent region are satisfied. We will — pgeenese. 
therefore write a second approximation \ J 
for f outside the laminar layer as ae J ‘ 
f = Bs + Bi(1—1) aa \ Pi 

+ Bi(1 — )?, (11.18) \. WA 








and satisfy the conditions at 7=1 and 
n=n*, where n* (see Fig. 6) represents Fic. 6. 

the effective thickness of the laminar 

sub-layer. Then, from (11,1), (11.2), (11.14), (11.16), (11.17), we have, when 7 =7*, 
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and when n=1, f=0f/dn=0. 
Hence, equating f and 0f/0n, as given by (11.18), to the above boundary conditions 
at 7 =7* and 7 =1, we obtain 


Be = Bi = 0, 


1 wsX 
, ws + bn*? {i+ sa (1+ ) n't 
Bi =—- - wane (11.19) 


2 (1 — n*) 








and 


1 3 wsaX 
1+ 0; + "% bn*? E + rin (1 + ~<) | 


a 
1 1 wsX 
+ ci | ont {1 aa Fz (: + ~~) nt - ws | — n*) = 0, (11.20) 
2 a 


which are the two equations from which to solve for Bz’ and n*. 
The distribution of the quantity (0f/0n) —w, in the turbulent layer is then given 
by the linear relationship 


of 
(= - a) - — ws — 2Bz(1 — »). (11.21) 
dn 


This is illustrated in Fig. 6, where the line A B represents the solution in the turbulent 
region (Eq. 11.21), and the curve OA the corresponding solution in the laminar layer, 
as given by Eq. (11.17). The chain curve OCB has been included for comparison, and 
is the solution developed in section 8 (b). The curves are approximately to scale, and 
refer to the same flow conditions. The effect of viscosity near the surface evidently 
has a large influence on the initial distribution of the shear stress, and hence upon 4;. 
For small values of 7 there is no reason to doubt the accuracy of the curve OA, but 
further from the surface the treatment of section 8 (b) is probably a closer approxima- 
tion than that afforded by the series (11.18). The true solution would therefore appear 
té be represented by a transition curve linking OA as 7-0 with OCB as n-1. 

One final point needs consideration. When 7=0, so also (0f/0n) —w;=q=0, and 
therefore the initial value of 00/8n becomes indeterminate. This may be circumvented 
by the following argument. From (11.16) 


Likewise, 








«wae ae 
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2 
Pg (*) “= R,X*n?. 
70 q1 


L (=) b " 
—j)= = — aX. 
70 On RiX? 


12. Solution of Eq. (5.9) for ’dH/ds. Since, by definition, 


=) 
We = Wel] > 
6 
we have, according to Eq. (9.3), 


ra a? (C _ *) (12.1) 
Wye = Wi B H? . ° 


Combining, now, Eqs. (4.6), (8.29), (8.53), (12.1) with Eq. (5.9), 
1-2 (8) fat afi ent Ee) 
_: iin Ao m ” B H? 


_— 
+f e ' re (12.2) 


0 A\? 
(1 + 1X ~) 
Ao 


where Js, a, 8 are found as described in section (10), and (0f/8n) —w,; is given by Eq. 
(11.17) in the range 0<<7*, and by Eq. (11.21) in the range n* <7 <1. 

Considering the functional nature of the terms, excluding H for the moment, on 
the right hand side of Eq. (12.2), it is evident that they are either functions of ws 
and R;, or of ws or R; separately. Hence we may write 


Hence, from Eq. (4.5), 














dH 
v ,% = ¢1(w;, R;, H). (12.3) 


Actually, from Eq. (9.4), H is also a function of ws and R;, so that (12.3) reduces to 


dH 
o— = $2(ws, R;). (12.4) 
ds 


For the numerical integration of Eqs. (5.9), (5.10) along the lines of refs. (1) and (2), 
however, it is desirable to retain the form (12.3). 

We now show how (12.2) may be evaluated in terms of the data discussed in sec- 
tion (6), namely, when the static pressure and the initial values of ® and H at each 
step are known. 

From Eggs. (9.3), (9.4) we have 








212 W. S. COLEMAN [Vol. V, No. 2 








vH vH 
$= = ; (12.5) 
aX A 
aXo — 
Ao 
or, 
H 
we = (wy)o ——————_- (12.6) 
A\é 
= 
Ao 
where 
v dp ae 1 dq 
(we)o = ae | oe 
pgiXo ds Xo gi ds 


Eq. (12.5), or (12.6), can be solved as follows. Given 8, Xo is determined from the 
law relating Xo and Ry, i.e. Falkner’s power law,® or.the logarithmic law of Squire 
and Young.® Now assume a value of ws; hence we obtain A/Apo (which is purely a func- 
tion of ws) by calculation according to section 8 (b), or, alternatively, from Fig. 3. 
X then follows. This enables a value of 6 to be found consistent with the assumed 
value of w;, and the corresponding value of R; is calculated. Finally, knowing both ws 
and R;, we obtain a from section (9), and, since His given, a second approximation 
for 6 (or ws) follows from Eggs. (12.5), (12.6). Hence, by trial and error, or by pre- 
established families of curves, we may solve for the relation between 6 and #. By vir- 
tue of Eqs. (12.5), (12.6), therefore, (12.2) is expressible in terms of 3, H and the pres- 
sure distribution, so that formally we have 

dH 
§ — = ¢$3|(wo)o, Ro, H]. (12.7) 
ds 

In conclusion it is worth noting that (wg)o is a multiple of Garner’s variable I’, as 
it is of von Doenhoff’s corresponding term. Thus, Garner takes 

Pilon (12.8) 


T= iS Ses 


q1 ds 
and assumes Falkner’s power law for the skin friction, viz. 
2 2, —1/6 
Xo = 0.006534Ry , (12.9) 


so that combining (12.8), (12.9) 





oe 1 dq 
Tr = 0.006534 —- — ’ 
Xo qi ds 
= — 0.006534(ws)o. 


Similarly, von Doenhoff’s term, which we will denote by I’, is 
wag 2 1 as 
Xo gi ds 


nis 2(we)o. 








th SOOT aT 
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Consequently, as discussed in section (6), both the theoretical Eq. (12.2) and the 
empirical equations for dH /ds of refs. (1) and‘(2) contain the same basic parameters. 

13. Comparison of Eq. (12.2) with the corresponding relations of von Doenhoff 
and Garner. For comparison we have taken, as examples, the case (a) when w;=0, 
(b) when w; = 30, the range of log Rs considered being in both cases 3.5 <log Rs <6.0. 
For (b) this range of Reynolds number allows a wide variation of flow conditions to 
be studied, from virtually the plane flow state up to separation, which is imminent 
when H exceeds 1.8. 























TABLE 1 
ws =0 
rot | | ‘ | odH/ds X10 
Pr a B xX H o/s | se | 
ai hal Qo/Y | | Ref. 1 | Ref. 2 | Theory 
3.5 | 2.845 | 20.21 | 0.0526 | 1.596 | 0.0937 | 2.471 | —0.996 | —2.731 | —0.047 
4.0 | 2.724 | 17.29 | 0.0464 | 1.421 | 0.0896 | 2.952 | —0.190 | —0.102 | +0.306 
4.5 | 2.682 | 16.10 | 0.0413} 1.330 | 0.0835 | 3.422 | —0.041 | +0.179 | 0.298 
5.0 | 2.667 | 15.57 | 0.0371 | 1.277 | 0.0776 | 3.890 | +0.007 | 0.202 | 0.271 
5.5 | 2.662 | 15.31 | 0.0338 | 1.242 | 0.0727 | 4.362 | 6.027/ 0.182] 0.190 
6.0 | 2.660 | 15.23 | 0.0309 | 1.215 | 0.0677 | 4.831 | 0.038 | 0.155 | 0.157 
| | A 
TABLE 2 
5=30 
' | | rm, edH/ds X10 
-0g a e | ¥ | 2 | o/s | “OB |——____ - 
q5/v | | | a0/Y | Ref.1 | Ref.2 | Theory 
3.5 | 7.770 | 103. a 0.0379 | 2.020 | 0.1460| 2.665| 35.6 | 72.7 | 62.0 
4.0 | 8.086 | 114.81 | 0.0334] 1.902 | 0.1421] 3.153] 20.7 | 35.1 | 40.9 
4.5 | 8.187 | 119.66| 0.0297] 1.770 | 0.1376 | 3.639 | 12.65| 15.4 | 24.6 
5.0 | 8.223 | 121.68 | 0.0267| 1.654 | 0.1328 | 4.123] 7.58| 7.58 | 14.9 
5.5 | 8.233 | 122.36| 0.0243) 1.566 | 0.1279] 4.607| 5.06] 4.34 | 9.23 
6.0 | 8.237 | 122.67| 0.0222| 1.495 | 0.1224] 5.088 3.60 | 2.62 | 5.68 








Tables 1 and 2 and Fig. 7 summarize the results, those of refs. (1) and (2) having 
been calculated from the formulae there given, which, for the present purpose, may 
be most conveniently expressed in the notation of the present paper as follows: 


Von Doenhoff and Tetervin' 


dH 
3 — = ef. 68(H-2-975) [2(ws)o — 2.035(H — 1.286) ], (13.1) 


ds 


with Xo given by Squire’s and Young’s formula 


= [5.890 logio (4.075Ry) |-?; 
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Garner? 


dH 2 
3 — = Xoe5(4-1- [(wy)o — 2.068(H — 1.4)], (13.2) 


ds 


where Xj is based on Falkner’s formula 


Xi = 0.006534R; 
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Fic. 7. Comparison of theoretical and empirical solutions for édH/ds. 


In estimating (ws)o, when ws; = 30, the skin friction law adopted in each case has been 
adhered to. Thus, for theoretical values of (wy)» Eq. (8.48) applies, whilst for the empir- 
ical formulae of von Doenhoff and Garner the particular skin friction equations given 
above have been used. This leads to slight variations of (w,)» under otherwise similar 
conditions, as will be seen from Fig. 7, where the appropriate value of (wy) is indi- 
cated against each point to facilitate comparison. 
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For plane flow (w;=0), both the theoretical and empirical results are in good agree- 
ment at large Reynolds numbers, but there is some discrepancy at lower values when 
H becomes abnormally large due to the very thick laminar sub-layer. Under these 
circumstances the theory is probably unreliable, so that some deviation is to be ex- 
pected. The condition that dH/ds +0, or, as is generally assumed, H =constant when 
the pressure gradient is zero or small, is, however, well substantiated by theory. 

Under conditions of a pronounced pressure gradient, tending to separation, agree- 
ment is less satisfactory, not only between the theoretical and semi-empirical solu- 
tions, but also between von Doenhoff’s and Garner’s results. 

Some comment on these discrepancies was made in the introduction. It was there 
pointed out that the theory is consistent with Garner’s relation in that both indicate 
ddH/ds to be a function of ws, Ry and H, whereas the von Doenhoff-Tetervin equation 
does not contain Ry as an independent. parameter. Consequently, in this respect, 
theoretical calculations might be expected to agree (as is the case) more nearly with 
Eq. (13.2) than with Ea. (13.1), the latter probably being less generally representa- 
tive on account of the above restriction. 

A further point of interest in regard to the empirical formulae is the value of H 
when dH /ds =w;=w3=0. From (13.1), (13.2) it will be noted that, according to von 
Doenhoff and Tetervin, H then has the value 1.286, whereas Garner, after investigat- 
ing the variation of H at transition in some detail, concludes that the value 1.4 is a 
better approximation. As already observed, 8dH/ds is very small within the normal 
range of H when a; is zero, so that any slight discrepancy to which the above quantity 
may be subject will introduce a large change in the value of H at which dH/ds is 
precisely zero. In the empirical approach dH/ds was obtained graphically, and there- 
fore the difference between the values of H under the afore-mentioned conditions, as 
given in refs. (1) and (2), is probably due to errors in the graphical method. By inter- 
polation of the theoretical results in Table 1, the corresponding value of H is found to 
be about 1.56. Allowing for the approximate nature of the theory, which may accord- 
ingly imply errors in the above H value of the same order as for the empirical 
formulae, it follows, nevertheless, that Garner’s figure of 1.4 is perhaps a better esti- 
mate than that of von Doenhoff and Tetervin. 

In conclusion, the general form of the equation for ddH/ds appears to be estab- 
lished, and a method of analyzing the growth of the turbulent boundary layer has 
been developed. In practice, however, it is laborious to use, though by extensive 
graphical treatment it is considered that the work could be reduced to reasonable 
proportions. Alternatively, if some degree of empiricism is acceptable, a reasonably 
reliable, simple and rapid means of reaching a solution is possible on the lines of refs. 
(1) and (2). 


APPENDIX 
According to Eq. (11.8) 


T i ‘é re) 0 
ShOMOR® 
eg? Rs On\ i dn\qi/| an\ qu 

The first term on the right hand side of (A) represents the true viscous shear stress, 


whilst the second term represents the Reynolds stress. Where the laminar sub-layer 
merges with the fully turbulent layer, the former stress becomes vanishingly small in 
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relation to the latter. Let n; denote the value of 7 at which the viscous stress is vir- 
tually negligible, and let the ratio of Reynolds stress to viscous stress at this point 


be m. Then 
0/q 0/4 m oq 
{OH9-220 . 
On\ qi/ | On\ qu Rs On\ qi 


Also near the surface * 


A = Ay, (C) 


2 





and approximately we may write 
a9/4 
—(—]) = RX’. (D) 
On \ qi 


Hence, by combining (B), (C), (D) and re-arranging, we have 





Vm , 
a ” 
so that if 40 of section 8 (c) is regarded as proportional to m, 
1 
nm~u™~ RX ; (F) 


By writing ~/m/A=constant, Eq. (E) becomes identical with von Kérm4n’s equa- 
tion” for the non-dimensional thickness of the laminar sub-layer. For the constant 
von K4rm4n finds the approximate value 30. Then if we take the condition of plane 
flow, when 

A= Ao = 0.392, 
we see that 


/m = 30 X 0.392 = 11.76. 


or 
m = 138.3. 


Thus the viscous stress is locally 1/138.3 of the Reynolds stress, i.e. it amounts to 
slightly less than 1% of it. 
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—NOTES— 
THE PRODUCT OF THREE REFLECTIONS* 


By H. S. M. COXETER (University of Toronto) 


. 

1. Introduction. This paper completes a trilogy, as it is a sequel to Synge’s “Re- 
flection in a corner formed by three plane mirrors” and Tuckerman’s “Multiple re- 
flections by plane mirrors.”! However, it is written in such a way as to be self-con- 
tained. 

The product of reflections in two intersecting planes is easily seen to be a rotation 
about the line of intersection of the planes through twice the angle between them. 
Hence the product of three reflections is, in general, a rotatory reflection: the product 
of a rotation and a reflection. We shall see (in §2) that, by a suitable displacement of 
the mirrors, without altering the final effect, we can arrange for the axis of the rota- 
tion to be perpendicular to the plane of the remaining reflection, so that the rotation 
and reflection are commutative. [An exception arises in the case of three “vertical” 
mirrors, all perpendicular to one plane. Then the product is a glide-reflection: the com- 
bination of a translation and a reflection. We shall discuss that case in Sec. 8. ] 

If the three mirrors have a common point, we represent them by the great circles 
which they cut out on the unit sphere around the point. Of the eight spherical tri- 
angles formed by the three great circles, just one is distinguished by having the three 
“fronts” of the mirrors all facing inwards. We shall call this the mirror triangle and 
denote it by ABC. [Similarly, if the three mirrors are all vertical, we represent them 
by the sides of the plane triangle ABC which they cut out on a horizontal plane. | 

Let AA;, BB,, CC; be the altitudes of triangle ABC; then A,B,C; is called the pedal 
triangle (or “orthic triangle”) of ABC. (See Fig. 1.) We shall find that the product of 
reflection in the three sides of an acute-angled triangle ABC, in any particular order, 
is a rotatory reflection [or glide-reflection | consisting‘of the reflection in one side of 
the pedal triangle combined with a rotation [or translation] along that side, of an 
amount equal to the perimeter of the pedal triangle. The three sides of the pedal tri- 
angle, with two senses along each, account for the six possible orders of the three 
reflections. We shall reconcile this with the results of Synge and Tuckerman, and 
finish with a brief statement concerning the analogous problem in m dimensions. 

Steiner remarked in 1837 that the pedal triangle of an acute-angled spherical tri- 
angle ABC is the triangle of minimum perimeter inscribed in ABC. The present treat- 
ment was suggested by Schwarz’s beautiful proof of the corresponding theorem for 


plane triangles.? 
For ‘simplicity, all the spherical triangles in the diagrams have been drawn as 


* Received Sept. 11, 1946. 

1 This Quarterly 4, 166-176 (1946), 5, 133-148 (1947). 

2 Steiner, Aufgaben und Lehrsdtze, Gesammelte Werke, vol. 2, Berlin, 1882, p. 45 (No. 7). Schwarz, 
Beweis des Satzes, dass unter allen einem spitzwinkligen Dreiecke eingeschriebenen Dreiecken das Dreieck 
der Hihenfusspunkte den kleinsten Umfang hat, Gesammelte Mathematische Abhandlungen, vol. 2, Berlin, 


1890, pp. 344-345. 
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plane triangles. This economy is justifted by the fact that the results remain valid in 
the case of three vertical mirrors. 
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The author wishes to express his gratitude to Professor Synge for many helpful 
suggestions in connection with this work. 

2. The sides of the pedal triangle. If a congruent transformation of the sphere 
leaves a great circle invariant (as a whole, but not necessarily point by point) it is 
either the reflection in the plane of the great circle, or a rotation about the perpendicu- 
lar line (which we call a rotation along the great circle), or the combination of both: 
a rotatory reflection. In the last case we shall call the great circle the line of action 
of the rotatory reflection. In general it is unique: no other great circle is invariant 
under the given transformation. The only exception is when the rotation is a half- 
turn, or rotation through 7; then the rotatory reflection reduces to the central inver- 
sion, which reverses all vectors and transforms each point of the sphere into its 
antipodal ‘point. 

Let A, B, C denote the reflections in the sides a, b, c of the mirror triangle ABC. 
The product Bc (first B, then C) is a rotation through —2A about the vertex A (or 
about the line to A from the center of the sphere). This can just as well be expressed 
as the product of reflections in any two great circles through A making this same 
angle A, say b’ and c’. Choosing b’ to be the altitude AAj, as in Fig. 2, we have 
ABC =AB’C’, where AB’ is a half-turn: the product of reflections in two perpendicular 
great circles through A). But this half-turn can just as well be expressed as the product 
of reflections in any two perpendicular great circles through Aj, say a’ and b’’. Choos- 
ing a’ to be A,Q perpendicular to c’, we have 4B’c’ =A'B"'c’, where a’ is perpendicu- 
lar to both b’’ and c’. Thus 4Bc =A’B’’c’: the product of the reflection A’ and the 
commutative rotation B’’c’. 

This shows that the line of action of the rotatory reflection ABC passes through Aj. 
Similarly that of CBA passes through C;. But ABc and CBA, being inverse transforma- 
tions, have the same line of action (in opposite senses). Therefore this line of action 
is the great circle A;,C;. Moreover, the sense of the rotation B’’c’ is from A; towards C,. 
Hence 


THEOREM 2.1. The six rotatory reflections ABC, CBA, BCA, ACB, CAB, BAC consist of 
reflections in the sides of the pedal triangle A,B,C, combined with rotations along those 
sides in the senses A,C,, C;A;, B,A;, AiB;, C:B;, B;C. 
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The amount of rotation, being twice the angle between b’’ and c’, is equal to 
2A,:Q, which is less than m@ for a sufficiently small triangle ABC, and consequently 
(by continuity) less than m for any triangle except a trirectangular one; for that is the 
only case where it can attain the value (which gives the central inversion). 

3. Two geometrical theorems. The rearrangement of reflections that led to the 
equation ABC =A’B’’c’ provides a kinematical proof for 


THEOREM 3.1. Jf A,B,C; is the pedal triangle of a spherical triangle ABC, the per- 
bendicular from A to A,C; makes an angle A with the altitude AA\. . 


Similarly the perpendicular from A to A,B, makes this same angle A with AAj. 
(See Fig. 3.) Thus AA; bisects the angle between the perpendiculars to A;C; and A;B;. 
Consequently it also bisects the angle between A;C; and A,B, themselves. Thus AA; 
is the internal bisector of ZB,A,C;, and BC is the external bisector. Hence 








Fic. 3 Fic. 4 


THEOREM 3.2. The orthocenter and vertices of a spherical triangle ABC are the in- 
center and excenters of its pedal triangle A,B,C,. 


This theorem can hardly be new, as its plane analogue is so familiar (see Fig. 1); 
but some of the standard textbooks fail to mention it. We have consideréd an acute- 
angled triangle ABC; but the effect of an obtuse C is merely to make C the incenter 
of A;B,C,; (so that the orthocenter of ABC is an excenter of A,;B,C;). 

The converse is, of course, obvious: If the external bisectors of the angles of any 
triangle A,B,C, form a triangle ABC, the internal bisectors are the altitudes of ABC. 

4. The perimeter of the pedal triangle. (For this part of the work we cannot allow 
any of the angles A, B, C to be obtuse.) By Theorem 3.2, Ais the center of an excircle 
of triangle A;B,C, (viz., the excircle beyond the side B,C,). The arc AiQ (Fig. 3), 
being one of the tangents from A, to this small circle, is equal to s:, the semiperimeter 
of triangle A,B,C;. As we saw at the end of §2, this arc is half the amount of the rota- 
tion B’’c’. Hence 


THEOREM 4.1. The amount of rotation involved in the product of reflections in the 
sides of an acute-angled spherical triangle ABC, in any order, is equal to the perimeter 
of the pedal triangle A,B,C,. 
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The right triangle AA,Q has side A\Q=s, and hypotenuse AA; =h,. Hence?’ 
sin s; = sin A, sin A 
cosS 2sinA/2 


=e (4.2) 
tan R tan R 





=2N=-—2 


where A is the area (or spherical excess) of triangle ABC, and R is the circumradius. 

The equation sin A,Q =2N remains valid when C is obtuse. But then (if A and B 
are still acute) the excircle of A,B,C; that has center A is the one beyond C,Aj, so the 
tangent A;Q is not s but s;—q, and the amount of rotation is 


2A,Q = 2(s1 - C1) = dj + by — ¢o, 


The following alternative proof of Theorem 4.1 is possibly of some interest. Be- 
ginning with the acute-angled triangle ABC, let us reflect in the side BC to obtain 
A’BC, then reflect the latter in its side CA’ to obtain A’B’C, and finally reflect 
A’B’C in A’B’ to obtain A’B’C’, as in Fig. 4. Let the respective pedal triangles be 
A;B,C,, A,B.Co, AoB:Cs, A;B;C3. Since BC is the external bisector of ZB, A)C,, the 
great circle C,A; contains B,; similarly it also contains C3 and Aj. The rotatory reflec- 
tion CBA, which transforms the triangles ABC and A,B,C; into A’B’C’ and A;B;C;,, is 
thus clearly exhibited as the reflection in the great circle C,A; combined with a rota- 


tion alung that great circle, of amount 
CiC; = A,A; = AiB, + B.C; + CsA; = ¢1 + a; + By = 251. 


On the other hand, when C is obtuse (while A and B are acute) the point B, is 
situated between C; and Aj, so we have 


CC; = AiA; = — B2A; + B.C; + CsA3s = —c1 +a,+b1 = 2(s1 — C1). 


5. The principle of transformation. Let S denote the sphere in any particular 
position. Let £ and F be two congruent transformations which change S into S¥ and S? 
(without moving the center). Then F changes S and S# into S¥ and S#¥. Since 
S#F =(S?)¥ ‘EF this means that F transforms the operation E into F—'zF. Whatever 
kind of operation E may be, F~'EF is the same kind; e.g., if E is a rotation through a 
certain angle about an axis /, F—'£F will be a rotation through the same angle about the 
transformed axis 1”. 

Now, since the reflection A transforms the rotatory reflection ABC into A~!-ABC-A 
= BCA, it follows that A reflects the line of action A,C; of the former into the line of 
action B,A; of the latter; so we have an alternative proof that BC is the external bi- 
sector of the angle B, A,C;. 

The same principle would have enabled us to foresee that the amount of rotation 
involved in ABC is the same as in BCA or CAB. (The remaining products CBA, ACB, BAC 
obviously involve the same amount of rotation, because they are the inverse opera- 


tions. ) 
6. The case of a right triangle. When two of the mirrors are perpendicular, the 


pedal triangle collapses. 





* For this use of the letters N and S, see Todhunter, Spherical trigonometry, London, 1914, pp. 29, 33. 
92, 109; or McClelland and Preston, A treatise on spherical trigonometry, London, 1886, Part I, p. 55, and 


Part II, p. 11. 
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If only one right angle occurs, one line of action is the altitude to the hypotenuse; 
to be precise, if C is a right angle (so that A; and B, coincide with C, and aB=BA), 
then CC, is the line of action of ABC=BACc (or, in the opposite sense, of CBA =CAB). 
By the principle of transformation, the only other line of action (viz., that of ACB 
and its inverse BCA) is the image of CC, by reflection in either of the two perpendicular 
mirrors. The amount of rotation is equal to the perimeter of the collapsed pedal tri- 
angle, viz., 2h,. 

If two right angles occur, there are only two inverse rotatory reflections, and the 
unique line of action is provided by that mirror which is perpendicular to both the 
others. But if all three angles are right angles, the product is the central inversion, 
and the line of action is indeterminate. 

To sum up: the number of distinct lines of action, regardless of sense, is equal to 
the number of oblique angles. 

7. Comparison with the work of Synge and Tuckerman. Synge represents the 
three mirrors by respectively perpendicular unit vectors Ny, Ne, Ns, and thence by 
points (having the same names) on the unit sphere. He shows that the product of the 
three reflections is a rotatory inversion, i.e., a rotation about a certain axis combined 
with the central inversion (p. 168). The amount of the rotation, denoted by @, is a. 
certain function of the vectors N;, Nz, Ns. The six possible permutations of the three 
mirrors lead to rotations through +8 about three axes, each combined with the cen- 
tral inversion. These three “optic axes” are represented by unit vectors A;, As, A3; 
and he shows very elegantly that the points Ni, Ne, N; are the midpoints of the sides 
of the spherical triangle A, AsA; (p. 170) and that 7+ is equal to the spherical excess 
of that triangle (p. 172). 

Since the central inversion may be regarded as a rotatory reflection involving rota- 
tion through 7, a rotatory inversion involving rotation through @ is the same as a rota- 
tory reflection involving rotation through +6 (or through +—@ in the opposite 
sense). The axis of the rotatory inversion joins the poles of the line of action of the 
rotatory reflection. 

Thus Synge’s N,N:N; and A; AeA; are the polar triangles of our ABC and A,B,C. 
His statement that N,, Ne, N; are the midpoints of the sides of A, A2A; is simply the 
dual of our statement that each side of triangle ABC makes equal angles with the two 
non-corresponding sides of triangle A;B,C;. (Theorem 3.2.) Similarly, dualizing the 
statement that AB is perpendicular to CC;, we find‘ that N; is distant 7/2 from the 
point of intersection (N,N:- A; A2). 

Since the angles of Synge’s triangle A,A2A; are supplementary to the sides of its 
polar triangle (our A;B,C,), the spherical excess of the former is equal to 


(w — aj) + (wm — 63) + (2 — 1) — oe = Qe — (1 +01; + 1) = 20 — 253. 


Thus 7 —@=2s , in agreement with our Theorem 4.1. 

Tuckerman obtains the same results again, representing the vectors by pure qua- 
ternions. His method shows very simply that +cos @/2 is equal to the scalar triple 
product of the vectors Ni, Nez, Ns (which Synge denotes by P). Consequently 


cos @ = 2P? — 1 = 8 sina sin (o — a) sin (¢ — az) sin (¢ — as) — 1, 


* Todhunter, op. cit., p. 114. 








222 NOTES [Vol. V, No. 2 


where a, a, a, o are the sides and semiperimeter of the spherical triangle NiN:Ns. 
Comparing Synge’s (3.10) with Tuckerman’s analogous Eqs. (13), we see that 


sin? 0/2 = k-? = cos? a; + cos? a2 + cos? az = 2 COS a COS a2 COS as. 
In terms of the mirror triangle ABC (with 2s;=a—6), this takes the form 
cos? s; = cos? A + cos? B + cos*C + 2 cos A cos B cosC. (7.1) 


Thus sin s;5=P=2N, in agreement with (4.2). 

8. Three mirrors forming a prism. When the common point of the three mirrors 
recedes to infinity, so that the mirrors are all perpendicular to one plane, we have 
essentially a two-dimensional problem: the product of reflections in the three sides of 
a plane triangle ABC. Work analogous to §2 gives AaBc=A’B’’c’, where A’ is the re- 
flection in a line a’ perpendicular to both b’’ and c’. We now have a glide-reflection: 
the product of the reflection A’ and the translation B’’c’ along the “line of action” a’. 
Theorems 2.1 and 4.1 remain valid with the word “rotation” changed to “translation.” 
The plane version of Theorem 3.1 is closely connected with the fact that the pedal 
triangle has angles r—2A, r—2B, r—2C. Finally, (4.2) reduces to the well known 
formula 


5s; = A/R 


for the semiperimeter of the pedal triangle of a plane triangle. 

9. The product of m reflections in m dimensions. Since the product of two reflec- 
tions is a rotation, the product of m is equivalent to [m/2] rotations (with an extra 
reflection if m is odd), and it is possible to arrange for the planes of the rotations to 
be all completely orthogonal, so that the rotations (and the extra reflection when it 
occurs) are commutative. Let the amounts of the rotations be 25:, 252, - - + , 25{m/2). 
We wish to express them in terms of the angles between pairs of mirrors. Let cj. 
denote the cosine of the internal angle between the jth and kth mirrors. Then® 


et2e1i gt2e2i ... (and —1 if m is odd) are the roots of the equation 
*t+1 — 2eyx — 2eig3x-- + — Winx 
— 262 x+1 — 2co3x +--+ — 2eomx | = O. 
_ 26m1 _ 2Cme — 2Cm3 “Le + 1 


When m=3, this reduces to 


(x+ 1)’ _ A(ces +- ca + cis + 2cescsiciz) x(x + 1) = 0; 


so we have 
/2 —1/2.2 


/ 2 2 
(2 cos s1) = (x +x )= 4(ces + 31 + C12 + 2C2s¢s1C12), 


in agreement with (7.1). 


———— ee 
* For a proof, see Coxeter, Lésung der Aufgabe 245, Jahresbericht der Deutschen Mathematiker- 


Vereinigung, 49, ¢-6 (1939) 
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AVERAGES OVER NORMAL MODES OF COUPLED 
OSCILLATORS WITH APPLICATION TO THEORY 
OF SPECIFIC HEATS* 


By ELLIOTT W. MONTROLL (University of Pittsburgh) 


1. Introduction. It is well known that the vibrational contribution to the specific 
heat at constant volume of a crystalline solid or of a gas composed of complex mole- 
cules can be calculated from the Einstein equation! 


C.= kd (hv;/2kT)?/sinh? (hv;/2kT), (1) 


where 
h=Planck’s constant 
k= Boltzmann’s constant 
T =Temperature in degrees absolute. 
| v;} =set of frequencies of normal modes of vibration. 
The summation is to be extended over all frequencies. Eq. (1) can also be written 
as the Stieltjes integral 


CL.=k f “hy/2kT)*4N(p)/sinh? (hv/2kT), (2) 


where N(v) is the number of normal modes with frequencies less than v; and vz is the 
largest frequency of a normal mode. : 

To apply Eqs. (1) or (2) one must either know all the frequencies of normal modes 
or their distribution function. In the case of complex molecules and crystalline solids 
composed of N (N=0O(102%) in the latter case) particles, the procurement of this 
knowledge becomes equivalent to the determination of the (3N —6) normal modes of 
N coupled oscillators (6 being the number of translational and rotational degrees of 
freedom). Mathematically this means finding the zeroes of a characteristic determi- 
nant of order (3N—6) or the roots of a (3N—6)th degree equation. In highly sym- 
metrical molecules and crystal lattices this characteristic determinant is factorable 
into several smaller ones; but to date N(v) has never been derived exactly for any 
molecular model of a real crystal. Numerical approximations have been obtained for 
N(v) but without any estimate of their errors.? 

By taking traces of powers of the matrix of the determinant mentioned above one 
can obtain the moments of the distribution function dN(v). It is the purpose of this 
note to show how, without a knowledge of N(v) or the frequencies of the individual 
normal modes, one can obtain by the method of mechanical quadratures averages 
over the normal modes 


* Received Dec. 4, 1946. 

1R. H. Fowler, Statistical mechanics, Comb. Univ. Press, 1936, pp. 86-103. 

2 E, W. Kellermann, Phil. Trans. Roy. Soc. 238, 513 (1940); M. Blackman, Proc. Roy. Soc. (A) 148, 
365-384 (1935); (A) 159, 417-431 (1937); P. C. Fine, Phys. Rev. 56, 355-359 (1939); M. Iona, Phys. Rev. 
60, 822 (1941); E. Montroll, J. Chem. Phys., 10, 218-229 (1942); 11, 481-495 (1943); E. Montroll and 
D. Peaslee, J. Chem. Phys., 12, 98-106 (1944). 








224 NOTES [Vol. V, No. 2 


F= 2 sod = f Hoare, (3) 


(of which (1) is a special case) from the traces of powers of the matrix of the charac- 
teristic determinant. An expression will be given for the upper bound of the error in 
the calculation of (3) from a finite number of moments of N(y). 

2. Résumé of theory of mechanical quadratures. In this section we shall quote 
without proof some theorems on mechanical quadratures that will be useful in our 
work. A detailed exposition of these theorems has been given by G. Szegé.?* 

Let a(x) be a non decreasing function with infinitely many points of increase in 
the interval (a, 6) and suppose the moments 


b 
is -{ x"da(x), n= 0,1,2,-:- (4) 
exist. 
If we orthogonalize the set of non negative powers of x, viz. 1,x,x?, - - + ,x" with 
respect to the function a(x) we obtain the set of polynomials po(z), pi(x),---, 
p»(x), +++ uniquely determined by the conditions 


a) p,(x) is a polynomial of precise degree m in which the coefficient of x" is positive 
b) the system {p,(x)} is orthonormal 


b 
f Pn(X) Pm(x)da(x) = bmn m,n=0,1,2,---. (5) 


Further, for n=1 
Ho M1 M2 ** Mn 
Fl Me Ms °° * Mn+1 
Pn(x) = (D,~1D,)—*/? Oa ae eee. ered. (6) 


| 
| Mn—1 Mn Mn+1 2 °° fen—1 
| 





1 x mt + 3 
Ho MA eo 


D, (7) 





| Mn Mn+1° °° Man | 


When a(x) has only a finite number, say m, of points of increase there exist only m 
linearly independent orthonormal polynomials.‘ 


Let x1<x%2.<%3< - ++ <x, denote the zeros of p,(x). There exist numbers 
Mi, Ae, * * * » An Such that by the Gauss-Jacobi mechanical quadrature formula 


b 
f F(x)de(x) = daF(21) + doF(22) + +++ + NaF (ae), (8) 


where F(x) is an arbitrary polynomial of degree 2n—1. The values of the numbers A 





3 Szegé, Orthogonal polynomials, Am. Math. Soc. Colloquium Ser., 23, Chapt. 2-3 (1939). 
4G, Szegi, loc. cit., 25-27. 
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of course depend on n. They are often called Christoffel numbers and are given by 
fhe following equivalent expressions 


id ee we Poe 
— J Fs (x)(x — =| doc(3), ed 





Rn+1 (-— 1) 
1 = — : ) (9b) 
kn Prn+i(Xy) Pn (x,) 
1 
My = {pola}? + {pila}? +--+ + {Pn(ar) }?. (9c) 
Here k, is the coefficient of highest power of x in p,(x) and from (6 
Re = (D,-1/D,)'". (10) 


Furthermore, all numbers X, are positive, and 
Ai + Ae +++ +An = a(b) — a(a). 


It has been shown by Markoff® that when F(x) is a continuous function with a 
continuous derivative of order 2n+1, (8) can be generalized to 


b 
[ F(x)da(x) = AF (a1) +++ + An (xn) + F(")(£)/(2n) !Rn (11) 


where a SED. Clearly, if (11) is approximated by (8) an upper bound to the error 
can be determined by finding the value of £ for which F@”(£) is a maximum and 
evaluating F(2”)(£)/(2n)!k?. 

We shall now show that the moments u, can be found in a vibration problem if 
a(x) is interpreted as the number of normal modes with frequencies less than v=xv , 
and therefore that (11) can be employed to find averages over all the nomal modes 
of vibration of the system. 

3. Theory of vibrations. When the kinetic and potential energies of a system of 
coupled elements are given by 


T = Weuhs +: - + uheds “ 
V = $(bugi + ++ + Bangn + 2br2gige + 2139193 + + + + + 2bn—1,nGn—19n)s (13) 


the frequencies v of normal modes are given by the roots of the characteristic equa- 








tion® 
(by;/432a?) — v? by2/4e?(aiade)'/? + + + Bin/40?(aya,)'/? 
bo1/4m?(a2a)'/? (boo/422a3) — v? ++ + bon /4a?(aea,)'/? 
| = 0. (14) 
| Bai/4m?(anaQy)'/? baa /4ar*(Gnde)!!? +++ (bnn/407a2) — v? 
If we let B be the matrix whose elements are b;;/4m*(a,a;)", and +m, tye, °°°, £¥ 


be the 2 roots of (14), and if 2y™ represents the sum of the mth powers of the roots, 


then we have 


5 A. A. Markoff, Differenzenrechnung, Leipzig, 1896, p. 29. 
6 E, T. Whittaker, Analytical dynamics, Camb. Univ. Press, 1927, p. 179. 
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zy = 2n, Sy2™t1 = 0, Zy?™ = 2 trace B*. (15) 


Further, if we interpret x in (4) as the ratio v/v, and a(x) as the number of values of 
v/v (where v is a root of (14)) which are less than x, it is clear from (4) and (15) that 


1 
Ho -{ da(x) = yv® = 2n, (16a) 
oul 
1 
Mom+1 -{ x?™*1da(x) = D(v/yz)?™t! = 0, (16b) 
= 
. —2m 
Mom -{ x?™da(x) = Z(v/v,)?" = 2v, ~—‘trace B”. (16c) 
at 


The averages of interest, (3), are f = Zf(v;), where the summation extends over all 
positive roots of the determinant (14). By defining a new function 


F(v) = f(r), if v20 (17a) 
and F(v) = f(—»v), if » <9, (17b) 
(3) can be written as 

F 1 . 
f = —f F(xv,)da(x), (18) 
27-1 


which is of the form of (11). Since the moments yp, can be determined from (16), 
p»(x) can be evaluated from (6) and f obtained by applying (11). 

4. Numerical example. To demonstrate the manner of application of (18) let us 
consider the steps required in the determination of the specific heat of tungsten. It 
is to be emphasized that the purpose of this example is merely to review the details 
of a calculation—too few moments are employed to obtain results that could be ex- 
pected to be in agreement with observed specific heats. 

Tungsten forms a body centered cubic lattice. Its equations of motion have been 
determined by Fine, and Montroll and Peaslee? who used the assumption that the 
Cauchy condition on the elastic constants is satisfied. The traces of the characteristic 
determinant were found in the latter paper and the moments of the frequency dis- 
tribution were determined. The moments, yu, there presented in table II differ some- 
what from those defined by eq. 4. Those of table II (when r=2/3), which we shall 
denote by u, are related to the u,’s of eq. 4, in a lattice of N atoms, by u,=6Nu, /v_z, 
so that the moments po - - - ws are 


3(6N) ’ 
(733/4608)(6N). 


uo = GN, Mi = Ws =u =--- = 0, He 
ws = (11/36)(6N), Me = (61/288)(6N), Ms 


Substituting these moments into (7) we obtain 
D)=6N, Di =}(6N)2, Ds = (1/36)(6N)?, 
Ds; = [65/2(36)*|(6N)*, Ds = [(65)(31)/16(36)5](6N)5. 
and from (6) 











ee 
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po(x) = 1/(6N)*", = pr(x) = 22x /(6N)*, — pax) = (3/2)2"/°(2a? — 1)/(6N)™?, 

p(x) = 2(2/65)!/2(18%* — 112x)/(6N)?/?, 
pa(x) = (1/4)(2/31)1/2(28824 — 306x? + 65)/(6N)!/2. 

The roots of p4(x) are 

x1 = — .8769, xe = — .5418, x3 = .5418, x, = .8769. 

Thus, from (9c) we See that 

Ar = Aq = (.2172)(6N), Ae = As = (.2828)(6N). 
Thus, choosing the function f(v) to be 
SQ) = f(xv1) = k(x6)*/sink? 6x, 
where 6=hv,/2kT we obtain from (1), (11), and (18 
C, ~ 3NR{ (.4344)(.87690)2/sinh? (.87690) + (.5656)(.54180)?/sinh? (.54186) }. 


Qualitatively our method is equivalent to replacing the entire frequency spectrum 
by a small number, say n, of specially chosen sharp frequencies. These frequencies 
and their weight factors are chosen so that the values obtained for all averages over 
polynomials of degree (2n—1), or less, are exact. In conclusion, it might be men- 
tioned that in general analogous methods can be used in evaluating averages over 
characteristic values of linear operators. 


THE SUBSONIC FLOW ABOUT A BODY OF REVOLUTION* 
By E. V. LAITONE (Cornell Aeronautical Laboratory) 
In cylindrical coordinates the Laplace differential] equation, which defines the ir- 
rotational incompressible fluid flow, becomes 


1 1 
bzz + brr + — or + — bee = 0, (1) 
Tr r 


where the last term vanishes when the flow has axial symmetry about the x axis. 
In this case a solution of Eq. (1) based on a source distribution f(x) per unit 
length along the x axis from x=0 to x=L is 


S is f(E)dé 
smn 








[@—or+epr’ 
1% fé)(x — &)dt 
= z= o a ’ 2 
i aia 4nJy [(x — £)?+ r?]3/2 (2) 
. ef f(é)dé 
v=9¢, = 





iJ, Teorey’ 


where v/u =(dr/dx)» satisfies the fixed boundary conditions given by the body shape. 


* Received Dec. 9, 1946. 
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The linearized differential equation for the velocity potential of a compressible 
fluid flow with axial symmetry is given by 


1 
B’bzr + Prr + — r = 0, (3) 
rT 


where for subsonic flow 
B= V/1 — M?>0. | 
The first method of Goldstein and Young! may be used to convert Eq. (2) intoa 
linear perturbation solution of Eq. (3) in the following manner 


a f(é)dé 
4 te ae re ’ 
. 4nB J o [(x — §)?+ (Br)? ]}/2 
1 rf = f(e(x — Hd 
A ee Ki ’ 
r md, [(x — &)? + (6r)?]3/2 


a r* S(&)dé 
aed a | [(x — 8)? + (6r)?]9/2 


Av dr Au 
(=)- G+ G2), 
thee | 0 dx] Un / 0 


Equation (4) can then provide a solution for a fixed given body shape for all Mach 


numbers less than unity (0<$31) as shown in ref. 1. 
If the substitution =x+rz is introduced into Eq. (4) and the Taylor Expansion 


is written as 








(4) 





fe) = f(a + Br) = Do eo" 7M), 


n= nN: 


then Eq. (4) becomes 
1 a (Br) wnt fied( g) (L—z)/Br 2"tldz 
F Worse ¢ 








Au = — — sieceemenneter 
4B p—0 n! z/Br (1 + 27)3/? 
© \n—1f(n) (L—z) /Br nd 
PO eS (5) 
4 nmo0 n! —2/6r (1 + 2?)3/2 
or, to the first order terms in Br, 
1 : 1 — 2x/L 
Au = — inB jae GE) + f’(x) [log (4*%/L)(1 — x/L) — 2 log Br/L — 2] 
+ f’"(x)-(L/2)(1 — 2x/L) + f(x) -(L2/12)[(1 — x/L)? + (x/L)?] + --- 
ao: oe g tds 1 
(™)(z)- ———|]+0 (6) 
+ fa) | a | tocnt 


1S. Goldstein and A. D. Young, The linear perturbation theory of compressible flow, with applications 
to wind tunnel interference, Brit., A.R.C., R. & M. 1909 (1943). 
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1 
av = — {2/(x)/Br — f"(2) [Ar log (6r/L)] + O(6r)}, @) 


T 


for 0<x<L. 
On the body surface to a first order approximation 


(Av/tUt.)o = (dr/dx)o = f(x)/2xBrou., (8) 
or 
f(x) = 2xBu,ro(dr/dx)o = Bu.S’, (9) 
where 
S= ate. 
Then 
Au 1 (S’ /1 — 22x/L Br x / x 
2 ) =—— —(- ~) - sv-(1 + log — — log 24/=(1 - =)) 
Use mw (2x \1— x/L L L L 
S§ 2x SM 2x 22? 
+ L(1- t= (1 - 244-0 f, (10) 
4 L 24 i, L? 


for any x other than 0 or L. 

Eqs. (8) and (10) provide the pressure distribution on (or along the streamlines 
r—0) a slender symmetrical body of revolution in either incompressible or subsonic 
potential flow. They also provide a first order Mach number correction for the sub- 
sonic potential flow on any body of revolution. For example, the surface pressure co- 
efficient (C,) at x =L/2 for a symmetrical body of revolution at any Mach number 
(M) less than unity would be given by 


Ce logx/1 — M? , 
( — ) wits . (11) 
Cpm—o max 1— (L*S /48S” + ut e's ) + log (r/L) 


Eq. (11) agrees with the expression obtained by Lees? for a slender prolate spheroid 


(S( =0 for n>2). 
The flow about any body of revolution is given to the first order by Eqs. (4) and 


(9) as 
2( a) mae f S'(O(x — Bat 
Ux (Oe 9 [(x — 24 (Br)2]3/2 


(=)-=f- S’(€)dé 
Ure 4n Jy [(x — §)? + (6r)?]*” 


Eq. (12) shows that at large transverse distances from a body of revolution 
f 








(12) 





- f “S'(8)-(x — Bab 
 deBiz? Jo ‘ 


Ue 
2 Lester Lees, A discussion of the application of the Prandtl-Glauert method to subsonic compressible 
flow over a slender body of revolution, NACA, TN 1127 (1946). 
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so that the subsonic flow wind tunnel wall correction for a body of revolution would 
vary as 1/(1—M*)8/2, 
Eq. (7) may be written 
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Therefore, for 0<x<L, 
f(x) = 2xu.Brr' [1 + (log Br/L)(3B8%rr!"/2 + Br*r)/2r + rr’! + 92) + O(r?) | 
= Bu,.S'(x) + O(r4logr), for B>O. (13) 


Eq. (13) shows that at least for incompressible flow (8 =1) Eq. (12) will provide the 
surface velocities accurate to the order r‘ log r=0. A corresponding statement cannot 
be made for subsonic flow however since the compressibility effects of the flow have 
been considered only to the first order in Eq. (3). 
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The pressure coefficient for compressible flow is given to the second order by 


is 


7 


— 2(Au/t,) — (Av/tn)? — (Au/te)*(1 — M*) + O(bx, $5), 


which in view of Eqs. (8) and (10) becomes 


Cp = — 2(Au/u.) — (r’)? + O(r! log r). | (14) 


Figure 1 shows a comparison of the exact incompressible potential flow pressure 
distribution on a prolate spheroid (ellipsoid of revolution) with that obtained from 
Eqs. (10), (12) and (14). Figure 2 shows the effect of Mach number on the maximum 
C, of various thickness ratio prolate spheroids as computed from Egs. (10) and (12). 
Figure 3 shows the effect of Mach number on the pressure distribution of a parabolic 
arc body of revolution as obtained from Eq. (14) in conjunction with Eqs. (10) and 


(32). 
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A CRITERION FOR STABLE CHARACTERISTIC EXPONENTS* 
By AUREL WINTNER (The Johns Hopkins University) 


1, Let f(t) be a real-valued, continuous, periodic function and let the unit of length 
on the ¢-axis be so chosen that the period of f(t) becomes 1. Then the linear differential 


equation 


a’ + f(t)x = 0 (1) 
has two solutions of the form 
x = eg(t), z= o'h(2), (2) 


where g(t) and h(¢) have the period, 1, of f(¢) and do not vanish identically. The num- 
bers +A occurring in (2), the characteristic exponents of (1), are determined mod 27i 
only. If \ is not a multiple of wi, it is clear that the two solutions (2) are linearly inde- 
pendent. If \ is a multiple of wi, it depends on the elementary divisors of certain linear 
substitutions whether (1) does or does not possess two linearly independent solutions 
of the form (2) (in the second case, the general solution of (1) contains a secular term, 
that is, a term containing the factor f). 

Let the differential equation (1) be called of stable type if every solution, x(t), re- 
mains bounded as t—>+ ~, and let the characteristic exponents +A of (1) be called 
of stable type if they are purely imaginary (including 0). Then it is clear that (1) 
cannot be of stable type unless the characteristic e »onents are, and that (1) must be 
of stable type if the characteristic exponents are of stable type and distinct. Needless to 
say, the two characteristic exponents, \ and —X, are considered as distinct if and only 
if they are distinct mod 277 (that is, if \ is not a multiple of 77). 

2. The usual way of calculating \ is supplied by Hill’s method of infinite determi- 
nants (the zeros of the resulting transcendental equation being all possible values, 


+ a, +A + 27, +rA + 4nt,--:, 


of the characteristic exponent). For certain purposes, a less elaborate procedure is 
more convenient. Quite a direct procedure happens to be contained in the classical 
proof (Cauchy-Fuchs-Floquet) for the mere existence of solutions of the form (2). In 
fact, this existence proof is based on the consideration of the linear substitutions of 
the “monodromy group.” These are binary substitutions of determinant 1, have a 
common pair of latent roots, and the logarithms of these roots are precisely the char- 
acteristic exponents. In other words, e** is the root of a reciprocal quadratic equation, 


say 
e+ — Jae*> + 1=0, (3) 


where @ is a constant uniquely determined by the periodic coefficient function of (1). 

Accordingly, the calculation of \ can be reduced to the calculation of a. But acan 
be calculated from f(t), by applying the method of successive approximations to those 
two solutions of (1) which are determined by the initial conditions 


* Received March 12, 1947. 
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x(0)=1, 2x(0)=0 and 2x(0) =0, x’(0) = 1. 


This procedure, which seems to go back to Liapounoff (cf. pp. 425-431 of vol. IV, 
part III (1902) of Forsyth’s Theory of Differential Equations), leads to the following 
representation of a: 


Put 
t 
F(t) = f (to) dt (4) 
0 
and 
‘ n—1 
F, = {F(1) — F(t) + F(t) } TT {F() — Fes) }, (5) 
k=1 
where »=1, 2,- +--+ (with the understanding that 


that is, if the product is vacuous); so that 
F,, = F,(h, et on, » ba). 


Then, if the constants a), @, - - + are defined by 


1} 4 tn-1 : 
an = - f at, f dty--- F(t, >> * , tn)dtn (6) 
2 0 0 “0 


(that is, if the integral of (5) over the n-dimensional polytope 
054,56185°': S451 (7) 
is denoted by 2a,), the series 
as21—a:+a—*:: +(— Ie t+-:- (8) 


is convergent and represents the coefficient in (3). This is the substance of the veri- 
fications which will be found loc. cit. 

This method does not seem to be generally known among applied mathematicians ; 
it is not mentioned in Strutt’s usually quite complete monograph. 

3. Since (6), hence (8), is real-valued, and since the roots of the quadratic equa- 
tion (3) are 

e® =a + (a* — 1), 

it is clear that the characteristic exponents +A are of stable type (that is, purely 
imaginary, including 0) if and only if 


—1 


a 


lA 
lA 


1, 
and that they are of stable type and distinct (mod 277) if and only if 


-—-1l<a<il. (9) 
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Hence, the stability of (1) depends on whether the sum of the numerical series (8) 
does or does not satisfy these conditions. 

Liapounoff has used this method in order to show that, if f(t) =f(t+1) is positive 
throughout, the characteristic exponents of (1) must be of stable type whenever the 
average value of f(t), that is, the constant term in the Fourier expansion of f(t), is 
sufficiently small; namely, less than 4 (cf. loc. cit.). Liapounoff’s proof can be inter- 
preted as amounting to an application of the second mean-value theorem in integral 
calculus. By using the first mean-value theorem instead of the second, two stability 
criteria of quite another type will be deduced in the sequel. Neither of them will 
assume that f(¢)20. The first of them is as follows: 

(i) Let f(t) be a real-valued, continuous function of period 1. Suppose that its average 


1 
kt -{ f(dt (10) 
0 


M = max | f()| (11) 


and tts absolute maximum 


satisfy the inequalities 
>) M*/(2n —1)!<yS1 (12) 


(which are compatible with min f(t) <0); for instance, that 


M*e™ < 6u <6. (13) 


Then the characteristic exponents of (1) are of stable type and distinct. 

That (13) is sufficient for (12), is seen by comparing coefficients. 

As to the possibility mentioned parenthetical] y after (12), it is sufficient to observe 
that, while (11) and (12) imply that 


u>O0, (14) 


(14) and (10) do not imply that f(¢) is non-negative throughout. 
4. According to (4), the first factor on the right of (5) is identical with the differ- 


ence 
J “fOdt — f : fiat, 


where 0 S/, $4 31. It follows therefore from (11) and (5) that 
n—1 
|F.| S$ MI]| F(t) — F(tess) |. 
k= 


On the other hand, from (4) and (11), 
| F(te) — F (tess) | S M(te — toss), 
where ¢¢41 St;, by (7). Hence, the inequality 


|F, | S M*D, 
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holds at every point of the polytype (7), if D.=D,(t, - - - , ts) is an abbreviation for 
the product 


n—1 
D, = [J (te — tes2), (16) 
kel 


the (non-negative) square-root of the discriminant of the polynomial | [?_,(z—t). 
The integral of the product (16) over the polytope (7) is known to have the value 
1/(2n—1)!. It follows therefore from (6) and (15) that 


| on | Ss 4M"/(2n — 1)!. (17) 


Hence, by the first of the inequalities (12), 
Dd | an | < du. (18) 
n=2 


According to (5) and the parenthetical remark which follows (5), the function F; 
is identical with the constant F(1). It follows therefore from (6) that a;=3F(1). In 
view of (4) and (10), this means that 


ay = dy. (19) 
Consequently, from (8) and (18), 
|a — (1 — 4u)| < 4x. 
In view of the second of the inequalities (12), the last formula line implies that 
0<a<il. (20) 
Since this, in turn, implies that (9) is satisfied, the proof is complete. 


The first of the assumptions (12) can be improved somewhat. In fact, it is clear 
from the proof of (15) that (15) can be refined to 


|F, | < M*D, — M(t; — tn)Dn, (21) 

where (t;—t»)D, is positive (within the polytope). 
5. Since only (9) was needed but (20) was deduced, it is natural to ask after a 
condition which, in contrast to (12), leads to the a-range complementary to (21), 


that is, to the range 
~1<as0. (22) 


Such a condition is contained in the following dual of (i): 
(ii) The pair of assumptions, (12), of (i) can be replaced by 


> M*/(Qn —1)!<4—ypS 2; (23) 
n=2 : 

for instance, (13) can be replaced by 

M*e™ < 24 — 6y S 12. (24) 


In fact, (17) and (19) were obtained without any hypothesis and are therefore 
applicable whether (12) is assumed or not. But (17), (19) and (8) imply that 
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|oa| s|1— 44] +4>0 M*/(2n — 1)1. 


Hence, by the first of the inequalities (23), 
Ja| <|1-—4u] +24. 


On the other hand, the second of the inequalities (23) can be written in the form 
up =2, which means that 


}1—4u| =4u-1. 


' 


This completes the proof, since the last two formula lines imply the inequality |a| <1, 
which is (9). 

Conclusion. If u, M are defined by (10), (11), then either (12) or (23) [and so, in 
particular, either (13) or (24)] is sufficient for stability. 

As an illustration, let 


f(t) = (a+ bcos 2nt)“', where 0<b <a; 


so that (1) becomes the equation known from the problem of frequency modulation. 
In this case, (10) and (11) reduce to 


pw = (a? — 5?) and M = (a — 5), 


and so the above inequalities supply explicit conditions for pairs (a, b) which are sure 
to be of stable type. Neediess to say, the resulting inequalities for a and b are just 
sufficient for stability. Incidentally, since f(t) is now positive, Liapounoff’s criterion, 


u <4, also is applicable. 


LOWER BUCKLING LOAD IN THE NON-LINEAR BUCKLING 
THEORY FOR THIN SHELLS* 


By HSUE-SHEN TSIEN (Massachusetts Institute of Technology) 


For thin shells the relation between the load P and the deflection ¢ beyond the 
classical buckling load is very often non-linear. For instance, when a uniform thin 
circular cylinder is loaded in the axial direction, the load P when plotted against the 
end-shortening € has the characteristic shown in Fig. 1. If the strain energy S and the 
total potential g=S—FPe are calculated, their behavior can be represented by the 
curves shown in Figs. 2 and 3. It can be demonstrated that the branches OC and AB 
corresponds to stable equilibrium configurations and the branch BC to unstable equi- 
librium configurations. The point B is then the point of transition from stable to un- 
stable equilibrium configurations. 

It was proposed by the author in a previous paper’ that the point A was the criti- 
cal point for buckling of the structure under external disturbances, using the S, € curve 
for “testing machine” loading and the g, P curve for “deadweight” loading. The load 
P for the unbuckled configuration of the shell corresponding to the point A was called 


* Received April 2, 1947. 
1H.S. Tsien, A theory for the buckling of thin shells, J. Aero. Sciences 9, 373-384 (1942). 
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the lower buckling load of the shell. The energy represented by the vertical distance 
from the point A to the curve BC is then the minimum external excitation required 


to cause the buckling at point A. 


However, if the external excitation is large, there is no reason why buckling can- 
not occur at the point B’ directly under the point B. The minimum external excitation 
required is then given by the energy represented by the distance B’B. This amount of 
energy is actually absorbed by the structure during buckling. Since the curve BA 
represents the final state of the structure after buckling, for buckling to happen be- 
tween B’ and A, energy is absorbed, and for buckling to happen between A and C, 
energy is released. But in any event, the lower limit of buckling load is definiteiy 


C CLASSICAL BUCKLING 








€ 


Fic. 1. 


given by the point B’, not the point A. Therefore 
the lower buckling load should be the load P cor- 
responding to the point B’. 

By referring to Figs. 11 and 13 of the afore- 
mentioned paper, and assuming a square wave 
pattern, we find the lower buckling stress ¢ of thin 
uniform cylindrical shells under axial load to be 
given by 


o = 0.42Et/R 
for testing machine loading and 


o = 0.19Et/R 


for deadweight loading. The. corresponding val- 
ues under the previously proposed criteria are 
o =0.46Et/R and ¢=0.298Et/R for the two cases. 
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CORRECTIONS TO MY PAPER 


ON AN EXTENSION OF THE VON KARMAN-TSIEN METHOD 
TO TWO-DIMENSIONAL SUBSONIC FLOW WITH 
CIRCULATION AROUND CLOSED PROFILES* 


QUARTERLY OF APPLIED MATHEMATICS, 4, 291-297 (1946) 
By C. C. LIN (Brown University) 
Page 292, line 3: for “ts now an a parallel footing” read “is now on a parallel footing.” 
Page 293, line 8: for “in the region exterior to” read “in the region Ry exterior to.” 


Page 293, line 9: for “such that Ry” read “such that.” 


* Received Jan. 20, 1947. 


CORRECTIONS TO MY PAPER 


THE ANALOGY BETWEEN MULTIPLY-CONNECTED 
SLICES AND SLABS* 


QUARTERLY OF APPLIED MATHEMATICS, 3, 279-290 (1946) 
By RAYMOND D. MINDLIN (Columbia University) 


Equations (2.4), (2.6), (6.2), (6.3): replace (1-+»1) by £, in the coefficient of V?T. 

Equation (2.6): replace a by ay. 

Equation (6.2) is obtained from (6.1) and (2.4) and not from (6.1) and (2.6). 

The reference to Biot’s analogies at the end of Sec. 4 should read “analogies. . . 
between gravity loading and boundary pressures and dislocations and between 
thermal loading and dislocations.” 

The integration by parts following Eq. (4.6) is incomplete with the result that 
the following corrections should be made: 

Equations (4.9), (4.10), (7.5), (7.15): multiply the last integral by 2. 

Equations (7.13), (7.20): replace [(1—v1)(1—v2) —1] by [(1—»1)(1 —v2) —2]. 


* Received April 22, 1947. 


BOOK REVIEWS 


Mathematical Theory of Elasticity. By I. S. Sokolnikoff with the collaboration of 
R. D. Specht. McGraw-Hill Book Co., New York and London, 1946. xi+373 pp. 
$4.50. 


There are five chapters in this book, three of them dealing with the fundamentals of the theory of 
elasticity, one dealing with the torsion and flexure of homogeneous beams, and the last chapter dealing 
with variational methods in elasticity. . 

In the chapter on analysis of strain we find a careful discussion of the properties of infinitesimal 
strain including Cesaro’s proof for the sufficiency of the six equations of compatibility. A brief section is 
added on finite strain, showing the difference between the Eulerian and Lagrangian approach and giving 
selected references to contempory work. To this reader it seems regrettable that no mention is made 
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of the early work by Kirchhoff, Boussinesq and the Cosserats on this subject and no indication given with 
regard to the nature of the progress which has been made since then. 

There follows a chapter on analysis of stress, containing all the essential facts on this subject, The 
third chapter is devoted to the subject of stress-strain relations. Among the topics discussed are the gen- 
eralized Hooke’s law together with a description of types of elastic symmetry; various forms of the com- 
plete system of differential equations of the linear theory for the isotropic body; the strain energy function 
including thermodynamic considerations often not included in elasticity texts; Kirchhoff’s uniqueness 
proof for the solutions of the linear theory; and finally Saint Venant’s principle. 

In the foregoing three chapters the author uses his own distinctive symbols for stresses and strains, 
and tensor notation to the extent of writing for instance the three differential equations for the stresses 
in the form r;;,;-+F;=0. No use is made of results of the theory of tensors and in order to study specific 
problems the author restates the fundamental equations of the theory in the notation which is customary 
in the American literature. 

Nearly one-half of the book is devoted to the theory of torsion and flexure, mainly on the basis of 
Saint Venant’s theory. The classical solutions of the torsion problem for the ellipse, rectangle, equilateral 
and isosceles right triangle are derived in the usual way. About twenty pages are devoted to an introduc- 
tion of complex-variable theory, as a preparation to Mushelisvili’s method for the boundary problem 
for the two-dimensional Laplace equation. As applications of this method one finds solutions of the torsion 
problem for the cardioid, for a loop of the lemniscate, for the section bounded by two circular arcs and 
for the inverse of an ellipse. There follows a discussion of the membrane analogy, a section giving refer- 
ences to experimental work, and a section on the torsion problem for multiply connected cross sections, 
with the example of Bredt’s formula for the thin walled section without partitions. 

In order to treat the Michell theory of torsion of shafts with variable circular cross section there is 
included a section on orthogonal curvilinear coordinates (without the use of tensor notation): As an ex- 
ample the solution for the conical shaft is given. 

After this there is presented the Bessel-function solution for the torsion of the circular cylinder twisted 
by load distributions other than that of the Saint Venant theory, and also the relation between the torsion 
problem for the nonisotropic beam and for the isotropic beam. 

There follows a treatment of the flexure problem with emphasis on the importance of leaving arbi- 
trary the coordinate system in the cross section instead of choosing the system of principal axes through 
the centroid of the section. The center of flexure is defined as the load point giving zero mean local twist. 
Here one misses mention of the possibility of alternate definitions of this point and of some of the recent 
literature on the subject. A discussion of the center of twist would have been welcome also. 

Explicit solutions of the flexure problem are presented for beams which have as cross sections an 
ellipse, a rectangle, an equilateral triangle, the section bounded by two concentric circles, and the cardioid. 

It is somewhat of an anticlimax to have the work of this chapter conclude with a treatment of the 
elementary (“technical”) theory of bending of beams with straight axis, including the solution of some 
very simple examples. 

In the last section we find references to special extensions of the theory of torsion and flexure. Some 
topics are omitted such as the work on torsion with variable twist of thin walled sections, or of the prob- 
lem of the effective width of I, T and box beam flanges, or of the effect of transverse shear deformation 
on the bending of beams. 

The last chapter of the book is devoted to a study of variational methods. Following the example of 
Trefftz in the “Handbuch der Physik” there are given careful statements of the theorem of minimum 
potential energy (assuming infinitesimal strains and Hooke’s law), and of the corresponding minimum 
theorem for the stresses. In somewhat unusual order the principle of virtual work is derived from the 
principle of minimum potential energy and stated (without qualification!) to be equivalent to the latter. 
There follows a proof of the reciprocity theorem; an introduction to the direct methods of the calculus of 
variations with the example of the twisted bar of rectangular cross-section; and a study of error estimates 
for the direct methods. 

The final section of the chapter on variational methods describes the method of finite differences, 
including some recent developments, as applied to the two-dimensional Laplace equation. 

Within the limits which the authors have set themselves—limits which to the reviewer appear to be 
somewhat narrower than the title of the book indicates—their work can be said to represent an interesting 


addition to the literature on the theory of elasticity. 
E, REISSNER 
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An index of mathematical tables. By A. Fletcher, J. C. P. Miller and L. Rosenhead, 
McGraw-Hill Book Company Inc., New York and Scientific Computing Service 
Ltd., London, 1946. viii+450 pp. $16.00. 


The need of an extensive index of available mathematical tables has long been felt by workers in 

applied mathematics. The usefulness of such an index is felt, not only in making possible the numerical 
solution of a mathematical problem by the use of existing tables but also in influencing the choice between 
various mathematical formulations of a problem to take advantage of such existing tables. Moreover, as 
pointed out in the preface, the index will be useful also to makers of mathematical tables, in showing 
what has been done, thus avoiding duplication and stimulating work in the direction of filling the main 
gaps. : 
The work is divided into two parts. Part I contains an index according to functions. Part II is a 
bibliography listed alphabetically by authors. An alphabetic index to Part I is also included. While Part I 
is the cardinal feature of the book, Part II gives a concise reference list of books and papers containing 
tables or dealing with such things as errors, in particular tables or the bibliography of tables. 

Part I indicates references, cross-references, scope of these references, and eventual errors on such 
material as—arithmetical tables, mathematical constants, tables of algebraic and elementary functions, 
tables of the higher functions, tables of integrals, tables of solution of transcendental equations, tables 
for numerical computations by differences, etc. 

In appraising the task undertaken by the authors, one should consider the widely scattered character 
of the information collected, the enormous quantity of material accumulated, its duplication and its 
eventual reliability and finally, the fact that the authors did not confine themselves to mere indexing but 
have indicated in a great many cases the errors involved and accomplished a selection by discarding ma- 
terial of mere historical value. 

The preface states—‘“I hope that this index will be a useful servant to the many scientists the world 
over who are working to build up a better civilization. To them this work is respectfully dedicated.” The 
international character of the accomplishment will not escape those scientific workers who have tried to 
obtain copies of material which is scarce or available only in a geographically remote location. 

It is essential to continue work along the lines of this index in keeping it up-to-date and also perhaps 
to facilitate speedy access to the material which might be urgently needed by the reader and might other- 
wise be difficult to obtain. Encouragement and help in accomplishing these further aims seem to fall 
particularly well within the responsibility of the newly created United Nations Educational Scientific 


and Cultural Organization. 
M. A. Biot 


Contributions to the study of oscillatory times series. By M. G. Kendall. Cambridge, 
at the University Press; New York, The Macmillan Company, 1946. 67 pp. $1.75. 


This work is an outline of the application of certain common schemes of statistical analysis to a few 
sequences of the type encountered in studies of economic time series. The viewpoint is that of the statis- 
tician. 
The book is divided into seven chapters. The first is concerned with the qualitative discussion of meth- 
ods to be employed, and with the construction of three time series from random sampling numbers ac- 
cording to linear autoregression equations of the second order. Chapters two through six are concerned 
with the analysis of these series by the correlogram, the periodogram and the variate difference schemes. 
Chapter seven contains the author’s summary, in which he concludes, as one might suspect, that schemes 
of analysis not based on considerations of the generating processes give rise to meaningless results. An 
appendix contains tables for harmonic analysis. 

As an account of inadequacies of the methods discussed, the paper will interest statisticians, but 
mathematicians and physicists who have been interested in the general time series problem in the broader 


sense will notice the lack of reference to a large body of the work done in this field. 
J. KRUMHANSL 























